Oscillations in delayed positive feedback systems

Positive feedback loops exist in many biological circuits important for organismal function. In this work, we investigate how temporal delay aﬀects the dynamics of two canonical positive feedback models. We consider models of a genetic toggle switch and a one-way switch with delay added to the feedback terms. We show that long-lasting transient oscillations exist in both models under general conditions and that the duration depends strongly on the magnitude of the delay and initial conditions. We then show the existence of long-lasting oscillations in specific biological examples: the Cdc2-Cyclin B/Wee1 system and a genetic regulatory network. Our results challenge fundamental assumptions underlying oscillatory behavior in biological systems. While generally delayed negative feedback systems are canonical in generating oscillations, we show that delayed positive feedback systems are a mechanism for generating oscillations as well.


Introduction
Positive feedback loops are a prominent network motif in many biophysical and biochemical circuits.They appear at multiple different scales ranging from gene regulatory networks through hormone secretion.Well-known examples of this motif include platelet activation during blood coagulation, [1][2][3] the Ferguson reflex during human child birth, 4,5 the Cdc2-Cyclin B/Wee1 network for cell cycle regulation, [6][7][8] the interaction between the central nervous system and hypothalamus-pituitary-adrenal axis, 9,10 and the Mos/MEK/p42 MAPK cascade. 6,113][14] Beyond bistability, models of positive feedback often do not exhibit more complicated dynamics such as oscillations.A very common motif for the generation of oscillations is delayed negative feedback, which is a control motif used in homeostatic systems such as the pancreas and its regulation of blood sugar. 15,16Models of delayed negative feedback are well-known to generate limit cycles emerging from a delay-dependent supercritical Hopf bifurcation.The effect of temporal delay on positive feedback systems has received much less attention.This article is a first step in delineating the effect of delay on positive feedback systems.
We study the effect of delays on positive feedback models by incorporating explicit delay into feedback terms of two wellknown models: a genetic toggle switch and a one-way switch. 7,12,17oth models describe two proteins that undergo production governed by a Hill function and linear degradation.In parameter regimes that facilitate bistability, we show that these models exhibit long-lasting transient oscillations before reaching a steady state.These findings suggest a mostly unexplored mechanism for generating oscillatory behavior in a biological system.Under certain conditions, the duration of these oscillations is so long as to be indistinguishable from true limit cycle oscillations.From the viewpoint of biology, the transient oscillations could last longer than the lifespan of the cell or organism that generates them.
Our findings challenge fundamental understanding of mechanisms that drive oscillations in biological and biochemical systems.Historically, when oscillations were observed in an experimental setting, it was assumed that they were most likely caused by delayed negative feedback.We argue against this assumption.
We begin by giving a preliminary overview of the positive feedback systems that we investigate.For completeness, we also describe how delays facilitate pulsatile dynamics in delayed negative feedback systems.Thereafter, we incorporate explicit delay in positive feedback systems and describe the means by which long-lasting but transient oscillations emerge.models with time delay in the feedback are analyzed later in the manuscript.
2.1.1The toggle switch.The toggle switch is a canonical system describing mutual inhibition between biochemical species such as proteins 12,14,17,18 (Fig. 1A).This system can be described by the ordinary differential equations (ODE) ‡ where x(t), y(t) are concentrations of transcription repressors and a 4 0 is the maximal protein production rate for both species (assumed here to be equal).The parameter n A Z measures the cooperativity between the repressors 13 and is positive for the mutual inhibition case.Because each molecule inhibits the production of its inhibitor, the feedback provides positive feedback of x onto x and y onto y.That is, x inhibits its inhibitor y (and vice versa), so the cumulative effect is positive feedback of x onto x and y onto y.For n Z 2 the system admits three real steady states (x*,y*): two stable nodes and one saddle point.The existence of two stable steady states classifies system (1) as bistable, a feature that can emerge in positive feedback systems. 7,12Fig. 1A-C depicts dynamics of this system.The two stable nodes in Fig. 1B (black circles) are characterized by the dominance of one protein over the other x Ã o 4 y Ã À Á and the saddle point (black triangle) describes coexistence (x* = y*).The basins of attraction of each stable steady state are separated by the stable manifold of the saddle point, which for this system is W S I = {(x,y)|y = x}. 13he Hill function § used for the feedback is a modeling choice but is not necessary for the results in this work, which hold for more general functional choices for the feedback.With a feedback function g(u) that describes inhibition, we require that g(u) Z 0 for all u, g should be monotonically decreasing, and g -0 as u -N.Similarly, if g describes activation, we require g(u) Z 0 for all u, g should be monotonically increasing, and g -A A R Z0 as u -N, for some constant A. For n 4 0, it is straightforward to see that g(u) = (1 + u n ) À1 satisfies the inhibition conditions.Similarly, for n o 0, g(u) = (1 + u n ) À1 satisfies the activation conditions.
One alternative popular function used for modeling biochemical systems is the Boltzmann function, Taking n o 0 in system (1) transforms the toggle switch into a one-way switch: a biomolecular feedback system pertinent in cellular differentiation 7 that exhibits mutual activation (Fig. 1D) as opposed to mutual inhibition.Under certain parameter sets, the one-way switch is also bistable, but in this case the two stable nodes (black circles in Fig. 1E) are characterized by both state variables being at their maximal levels or completely extinct. 7,12,17The saddle point (black triangle in Fig. 1E) is characterized by both state variables coexisting somewhere in the middle.The stable manifold, as in the toggle switch system, runs through the saddle point and separates the basins of attraction of the stable nodes.Linear stability analysis around the saddle point yields, locally, a stable manifold equation W S A = {(x,y)|y = K À x}, with K E 1.36 when a = 2.15, n = À2.

Delayed negative feedback
A canonical example of a biomolecular feedback system is a protein inhibiting its own production after some temporal delay, t 4 0 (see Fig. 2A).Mathematically, this system can be represented by the delay differential equation (DDE) where a is the maximal protein production rate and n 4 0 characterizes the sensitivity of the negative feedback.Initial data for DDEs are prescribed by a history function x(t) = h(t) for t A [Àt,0].Although this system is comprised of a single equation, it contains limit cycles in its phase space that are generated by the delay.0][21] Furthermore, the frequency of the delay-induced oscillation decreases as t is increased and the amplitude scales as ffiffi ffi t p near the bifurcation point.Linear stability analysis of eqn (2) about the steady state x = x* yields the following conditions for a Hopf bifurcation: where o is the frequency of the resulting oscillation.Clearly, the Hopf conditions cannot be satisfied in the absence of delay (t = 0).Observe that taking t = 0 in eqn (2) renders it a gradient system, which cannot have oscillations.
For t 4 0, however, there exist critical values (o c ,t c ) satisfying the bifurcation conditions.When t is increased past t c , a pair of complex conjugate eigenvalues crosses the imaginary axis.Although this is not sufficient to guarantee the emergence of a stable periodic solution via a supercritical Hopf bifurcation for t 4 t c , the existence of stable oscillations beyond the Hopf bifurcation point can be verified numerically (see Fig. 2B and C).
We next discuss how the network motifs and dynamics presented in this section produce pulsatility in delayed positive feedback systems.

Delayed positive feedback
To investigate the role delay plays in producing oscillations in the toggle switch and one-way switch systems, we implement explicit delay, t, into both feedback terms of system (1).This yields the set of DDEs where all parameters are defined as in system (1).The cooperativity parameter, n A Z can be adjusted to model delayed mutual inhibition or delayed mutual activation by letting n 4 0 or n o 0, respectively.Importantly, the set of possible steady states for system (3) is identical to system (1).In all simulations, we use constant history functions for initial data.
In this section, we discuss dynamics of system (3) and how they change with t and different history functions.In particular, we discuss how the presence of delay facilitates onset of pulsatile dynamics.DDEs are infinite-dimensional structures, so illustrating dynamics in their true phase space is not possible.We therefore present results by projecting the phase portrait into a two-dimensional subspace of the phase space.

Delayed mutual inhibition
We set n 4 0 in eqn (3) to model delayed mutual inhibition.To see that oscillations are indeed present in the phase space of this delayed positive feedback system, observe that taking system (3) with initial data satisfying x(t) = y(t) = a A R for t A [Àt,0] constrains the dynamics to the stable manifold of the saddle point, W S I , which is the seperatrix partitioning the basins of attraction of the two stable steady states in the system without delay (eqn ( 1)).The resulting differential equation is exactly eqn (2), which is a model for delayed negative feedback and exhibits stable oscillations for sufficiently large values of t.Thus, eqn (3) admits pulsatile dynamics along W S I in the xyplane.Although this is a limiting circumstance, it illustrates that limit cycles are indeed possible in delayed positive feedback systems.In more complicated mutual inhibition systems, choosing history functions precisely on an invariant manifold is very difficult.For that reason, we investigate dynamics of (3) starting near, but not on, W S I .Numerical solutions suggest that limit cycles do not arise with initial data off W S I .However, the system exhibits transient but long-lasting oscillations that eventually contract to a steady state.Fig. 3B shows two sample trajectories of system (3) starting on opposite sides of W S I .Superimposed on top are the trajectories of the ODE system (1) with the same initial conditions.The ODE solutions (black) approach their respected steady states without oscillatory dynamics.In contrast, the DDE system trajectories (red and blue) oscillate for many cycles before eventually reaching their respective steady states.
Because the trajectory initially oscillates near the stable manifold of the ODE system, we refer to the long lasting oscillations as ''remnants of the delayed negative feedback oscillator'', which is present exactly on W S I .Numerical time series show that the period, T, of the transient oscillations is related to the delay by T E 2t.Thus, during the period of transient oscillations, , where N is the oscillation number.The time series in Fig. 3C shows oscillations in x graphed as a function of oscillation number.From a biological viewpoint, it is possible that the transient oscillations last longer than the lifespan of the cell or organism.In such a case, there is little distinction between ''long-lasting transient'' oscillations and ''stable limit cycle'' oscillations.There appear to be two main factors that contribute to the duration of the transient oscillations in the delayed mutual inhibition system: (1) the magnitude of the delay, t, and (2) the distance of the initial conditions from W S I .To quantify how these factors affect oscillation duration, we simulated the mutual inhibition system (3) and recorded the number of cycles in the x variable.The results are given in Fig. 4A.The number of oscillations increases super-linearly with t.
We next numerically investigated the delayed mutual inhibition system under varying initial conditions.To see how

Delayed mutual activation
We set n o 0 in system (3) to model delayed mutual activation (Fig. 3D).Similar to the previous model, this system displays bistability and sample trajectories on either side of W S A are shown in Fig. 3E.As in the delayed mutual inhibition case, we superimposed the trajectories from the ODE system (1) with identical parameters and initial conditions.The results show that the trajectories without delay (black) approach the stable steady states without oscillation, but the trajectories with delay (color) oscillate roughly parallel to W S A before eventually approaching a steady state.A sample time series for the state variable x is given in Fig. 3F and exhibits long-lasting oscillations before reaching a steady state.
As in the previous model, we simulated the delayed mutual activation model for varying delay values t to explore how oscillation duration changed.Results are summarized in Fig. 5A.For small t the system displays a small number of transient oscillations.However, increasing t to sufficiently large value leads  to a rapid rise in the number of oscillations, which continues to increase for larger delays.To understand the origin of this rapid increase in oscillation number, it is helpful to view the dynamics in the xy-plane (Fig. 6A and B).The trajectory in Fig. 6A corresponds to the parameters and delay value of the red point in Fig. 5A.With these parameters, the trajectory oscillates as it moves towards the upper steady state.However, when the delay is increased from 4 to 6, as in Fig. 6B, the trajectory first moves towards the upper steady state before changing direction to the lower steady state.This flow reversal generates many more oscillations, and there appears to be a critical delay value at which this dynamic behavior occurs.
In addition to varying delay values, we simulated the delayed mutual activation model with varying initial conditions.As with the previous model, we extended a normal line from W S A and sampled initial condition pairs along that line.Next, we simulated the system for each initial condition and recorded oscillation duration.The results for the number of oscillations as a function of the distance (D) to W S A are given in Fig. 5B.The plot suggests that, initially, as you get farther away from W S A , the number of oscillations increases monotonically.However, when D is sufficiently large there is a sharp downwards jump and the number of oscillations monotonically decreases as D is increased further.Looking again in the xy-plane (Fig. 6C and D) we plot two trajectories with identical parameter values, but starting at different distances from W S A .In one case (Fig. 6C), the trajectory begins to travel toward the upper steady state, but then turns around and crosses over W S A to eventually end at the origin.Initial conditions satisfying D o 0.07 in Fig. 5B show similar dynamics in phase space to that of Fig. 6C.Taking an initial condition a small distance farther away (green point) yields a trajectory that oscillates directly towards the upper steady state.Similarly, all starting points satisfying D 4 0.09 in Fig. 5B show similar dynamics in phase space to that of Fig. 6D.

Asymmetric networks
Thus far, we assumed that the dynamics of both of the interacting species are the same, described by eqn (3) with identical parameter values for x and y.We now explore the effects of breaking the symmetry, using the set of DDEs where the parameters values are not necessarily equal between species.Fig. 7 shows that long lasting oscillations can occur when the time delays differ between species, or the Hill coefficients differ, or both.Therefore, the phenomenon is not the result of symmetry between the interacting species.

Single equation example
The long-lasting oscillations in Fig. 3 and 7 come from models which describe the evolution of two species.However, we found that long lasting transient oscillations can occur in models of a single species.Consider the following DDE where a is the maximal production rate and n o 0 is the Hill coefficient.This equation describes a species which positively auto-regulates itself after some delay t.This system is capable of producing long lasting oscillations over time (Fig. 8B), similar to those produced by the model with two interacting species (Fig. 3F).

Biological examples
Here we present two, more complicated biological models of delayed mutual inhibition and delayed mutual activation.We show they are capable of displaying similar dynamics to system (3).

The Cdc2-Cyclin B/Wee1 system
Cdc2 is one member of a class of enzymes coded by celldivision-cycle (cdc) genes that is required to move through the cell division cycle. 22In preparation for progressing into the M phase of the cycle, cyclin B is synthesized and binds with free Cdc2 to form an inactive M-phase promoting factor (MPF).Only when this MPF is phosphorylated on threonine-167 by a cyclin-dependent activating kinase (CAK) is it in an active state.
Once the concentration of active MPF is sufficiently high, the cell can progress to the M phase. 23Wee1 is a protein that helps mitigate this process by phosphorylating the MPF at the tyrosine-15 site, rendering it inactive.Conversely, the MPF phosphorylates Wee1, rendering Wee1 inactive.This is a clear example of mutual inhibition, and mathematical models have shown this system to be bistable (high MPF, low Wee1 or high Wee1, low MPF). 7,8,24When Wee1 is high, the cell remains in interphase.
We extend the equations in a previous model for this system 6 to include explicit delays in the relevant feedback terms.The equations to model the interaction between the active Cdc2-Cyclin B and active Wee1 molecules are given by the dimensionless DDEs where x denotes the active Cdc2-Cyclin B and y represents the active Wee1.With the parameter values in Fig. 9, the system is bistable. 6In the system without delay, the basins of attraction of the two attractors are separated by the stable manifold of a saddle point.The exact equation of the stable manifold is not known for this system.As a proxy for it, we use a linear approximation to this manifold at the saddle point and denote it by W S Cdc2 .Fig. 9B displays a phase portrait of the system projected onto the xy-plane with sample trajectories (color) starting with initial conditions on either side of W S Cdc2 .The black trajectories are the solutions to the system without delay  starting at the same initial points.These trajectories exhibit similar behavior to those in Fig. 3B, demonstrating that it is possible to achieve transient oscillations in more complex delayed mutual inhibition models for biological systems.A sample scaled time series showing how active Cdc2-Cyclin B evolves is shown in Fig. 9C, exhibiting over 30 oscillations before reaching equilibrium, similar to system (3).

Genetic regulatory networks
6][27][28][29][30] Here we consider a simple delayed mutual activation model between a gene and protein shown schematically in Fig. 10A.A typical non-dimensional model for such a GRN, but incorporating a time delay in the feedback terms, is where G represents the fraction of the gene in its activated state and P represents the concentration of the gene product protein.
The parameters r, s, l, d, k, n represent the activation rate of the gene, inactivation rate of the gene, the production rate of the protein, the degradation rate of the protein, the half-max protein response, and the cooperativity, respectively.Under certain parameter sets, the system is bistable as shown in the phase plane in Fig. 10B.As in earlier examples, the dashed line is a linear approximation to the stable manifold (for the system without delay) of the saddle point, referred to as W S GRN .The trajectories of the system with time delay have a behavior that is similar to trajectories in Fig. 3E of the delayed one-way switch model.This illustrates that very long transient oscillations can occur in delayed mutual activation GRN models.A sample time course of protein concentration is shown in Fig. 10C.We observed long lasting oscillations in the protein concentration up until the trajectory converged to steady state.
The major implication of this result for GRNs is that the mechanism underlying protein production at the genetic level may manifest as a delayed positive feedback system.Historically, oscillatory protein production has been interpreted as a result of delayed negative feedback.With the sheer number of oscillatory genes in biology, our findings suggest the likelyhood that delayed positive feedback is an underlying mechanism in many cases.
As an example, gene expression of over 1000 proteins during development in C. elegans have been identified as undergoing oscillatory dynamics.Cuticular collagen genes, hedgehog receptors, and metallopeptidase, amongst many others, have been observed to oscillate. 31Since only a few cycles are typically observed, it is unclear whether these are true sustained limit cycle oscillations or if they are simply long-lasting.In the latter case, delayed positive feedback could be the driving mechanism.

Discussion
We presented two canonical examples for positive feedback, the toggle switch and the one-way switch, and explored their dynamics when adding temporal delay to the feedback terms.We showed that they are capable of generating long-lasting oscillatory dynamics under certain conditions.We verified our results using a two-dimensional projected phase space for comparison against the model dynamics without a time delay.Finally, we showed that qualitatively similar oscillatory dynamics occur when applying temporal delay to previously studied biological models of the toggle switch and one-way switch.Because of the large number of cycles produced, the transient oscillations in the delayed positive feedback systems may be indistinguishable from true limit cycle oscillations.Although our models used Hill function formulations for the feedback, we found similar long-lasting oscillations when either Boltzmann or hyperbolic tangent formulations were used (not shown).
We found that two main factors contribute to the duration of these transient oscillations: (1) the magnitude of the delay (t) and ( 2) the Euclidean distance (D) of the initial conditions from  W S I or W S A .The delayed mutual inhibition system showed that the oscillation duration increases super-linearly with delay and decreases super-linearly with D. The delayed mutual activation model showed oscillation duration increases with the delay value, but a sharp jump occurs at some critical t value.We showed in Fig. 6A and B the trajectories with t values on either side of the jump approach different steady states.We believe this may be due to a shift in the DDE seperatrix when changing the delay value t.This would explain why the trajectories approach different steady states even though the initial conditions are the same.The oscillation duration as a function of D shows initial increase in oscillation duration with a sharp decrease in oscillation duration at some critical distance.Similarly, we hypothesize that the critical D value in Fig. 5B could mark the location of the DDE seperatrix, and that on the true seperatrix a limit cycle exists.Thus, as initial data is taken further away from the critical D value, we observe fewer oscillations.We hope to explore these phenomena in the future to get a better idea of the shape of the true DDE seperatrix and determine the mechanisms responsible for these interesting transitions in the dynamics.
We observed that within parameter regimes which enable bistability, the long lasting oscillations in models of delayed positive feedback oscillated ''orthogonally'' to their equilibria.For example, with delayed mutual inhibition the equilibria are characterized by either x c y or x { y.However, x E y during most of the transient oscillation.In the delayed mutual activation case, the opposite occurs.The equilibria are characterized by both species completely extinct or maximally expressed.However, in the oscillation, the trajectories flip between states where x c y and x { y, which is reminiscent of mutual inhibition.This is significant, especially to experimentalists, because it gives an inference strategy for determining which mechanism produced the oscillations (if they occur due to delayed positive feedback).In-phase oscillations indicate delayed mutual inhibition, while anti-phase oscillations indicate delayed mutual activation.
2][33] Models for these behaviors would most likely contain a delayed negative feedback motif or a combination of negative and positive feedback, as these two systems are canonical generators of pulsatile dynamics.However, our results argue that the mechanism behind oscillations in such systems can potentially be explained by delayed positive feedback alone, as long as there is delay in the feedback and as long as the system has strong non-linearity.
Finally, although the examples shown here were based on biochemical or biological examples, the concept of oscillations driven by delayed positive feedback is applicable in other areas as well.Any mechanism that exhibits bistability must have an inherent positive feedback loop; thus, any system exhibiting bistability and strong non-linearity has the potential for long-lasting oscillatory dynamics.Examples include optical bistability 34 in condensed matter physics, physical biogeochemistry of sediments and oceans as a means to understand anoxic marine systems, 35 and astrochemical bistability arising from autocalatyic oxygen chemistry. 36

a b þ 1 ,
where a,b are constants.For b 4 0, this equation satisfies the inhibition conditions.Similarly, for b o 0, the equation satisfies the activation conditions.Another example is the hyperbolic tangent function, gðuÞ ¼ tanh x À where a, b are constants.When b o 0 (b 4 0), the equation satisfies the inhibition (activation) conditions.Using such functions would not alter the qualitative results of our work.2.1.2The one-way switch.

Fig. 1
Fig. 1 Reaction network diagram, phase space, and time series for system (1) when a = 10, n = 2 (A)-(C) and a = 2.15, n = À2 (D)-(F).Leftmost column: schematics of biomolecular feedback structures studied.Middle column: phase portrait of eqn (1) with sample trajectories (color).The dashed lines are the separatrices W S I (top) and W S A (bottom).Right column: sample time series for the trajectories shown in the middle panels.

Fig. 2
Fig. 2 Delayed negative feedback.(A) Schematic of the delayed negative feedback biomolecular system.(B) Oscillations emerge in system (2) as the delay is increased from t = 1 (blue curve) to t = 3 (orange curve).(C) Bifurcation diagram for eqn (2) when varying the delay t.Solid and dashed black lines correspond to stable and unstable stationary steady states, respectively.The red curve shows the minimum and maximum values of x on the periodic branch, and the green circle is a supercritical Hopf bifurcation point with t c E 1.7 and a = 10, n = 2.

3
Dynamics of eqn (3) for a = 10, n = 2, t = 6 (A)-(C) and a = 2.15, n = À2, t = 8 (D)-(F).Left column: schematics of the biomolecular feedback system.Middle column: dynamics in the xy-plane with actual or approximate separatrices W S I (top) and W S A (bottom) superimposed as dashed lines and sample trajectories from eqn (3) plotted (color).Black trajectories are solutions of eqn (1) for the same initial conditions.Right column: sample time series.The abscissa shows the oscillation number, using the fact that the period T E 2t to scale time.

Fig. 4 (
Fig. 4 (A) Number of oscillations as a function of t for the delayed mutual inhibition model.Initial data are x(t) = 4 and y(t) = 3 for t A [Àt,0] in each case.(B) Number of oscillations as a function of the distance from W S I .The delay t = 6 is each case.For both panels, a = 10 and n = 2.

5
(A) Number of oscillations as a function of delay value t for the delayed mutual activation model.Initial data are x(t) = 1.2 and y(t) = 0.5 for t A [Àt,0] in each case.(B) Number of oscillations as a function of the distance to W S A .The delay value is t = 6 for each case.For both panels, a = 2.15 and n = À2.

Fig. 6
Fig. 6 Dynamics of the mutual activation model depend critically on t and D. (A) and (B) Two trajectories with the same starting points (black stars), but with two different delays, t = 4 (A) and t = 6 (B) corresponding to the red and green points in Fig. 5A.(C) and (D) Two trajectories with the same delays (t = 6) but starting at two different distances from W S A , D = 0.07 (C) and D = 0.09 (D), corresponding to the red or green points in Fig. 5B.

Fig. 8
Fig. 8 Delayed positive autoregulation of a single species.(A) Schematic of biomolecular feedback structure.(B) Long lasting transient oscillations occur over time in system (5).Parameter values were a = 2.15, n = À2, t = 14.

Fig. 10 (
Fig. 10 (A) Schematic of the delayed gene/protein system.(B) Phase space with seperatix W S GRN superimposed (dashed line).Trajectories from eqn (7) plotted (color) along with those starting at the same location but without delay (black).(C) Sample time series showing extensive oscillations before reaching steady state.The parameter values are r = 0.1, s = 2, l = 9, d = 1, k = 2, n = 3, and t = 10.