# Pitchfork and Hopf bifurcation threshold in stochastic equations with delayed feedback

###### Abstract

The bifurcation diagram of a model nonlinear Langevin equation with delayed feedback is obtained numerically. We observe both direct and oscillatory bifurcations in different ranges of model parameters. Below threshold, the stationary distribution function is a delta function at the trivial state . Above threshold, at small , with at threshold, and monotonously increasing with the value of the control parameter above threshold. Unlike the case without delayed feedback, the bifurcation threshold is shifted by fluctuations by an amount that scales linearly with the noise intensity. With numerical information about time delayed correlations, we derive an analytic expression for which is in good agreement with the numerical results.

###### pacs:

02.30.Ks, 05.10.Gg, 05.70.Ln, 87.16.Yc, 87.18.CfWe obtain by numerical means the complete bifurcation diagram of a generic nonlinear and non Markovian Langevin equation that incorporates the effect of delayed feedback. Both pitchfork and Hopf bifurcation thresholds are observed and studied, and the results contrasted with two related limits: The deterministic limit of a differential delay equation, and the stochastic bifurcation of the same model without delay.

The study of differential delay equations re:driver77 () is an important topic in applied mathematics, with widespread applications in Physics (lasers, liquid crystals), control systems in Physiology (neural and cardiac tissue activity) re:mackey77 (); re:glass88 (), and Economy (agricultural commodity prices). Recent interest has arisen in the mathematical modeling of cellular function at the molecular level, especially in transcriptional gene regulation re:hasty01 (). Feedback regulation is a common motif in complex cellular networks, with delays arising from the complexity of the underlying network, or from the wide disparity in time scales of the many chemical processes involved in regulation re:glass88 (). For example, DNA is transcribed at a rate of 10 to 100 nucleotides per second, and it may take a delay of the order of minutes before the transcription factor appears as a finished product in the cell and is available for regulation. Significant delays can also be attributable to the time required for the diffusion of proteins through membranes, so that, for example, the auto regulated feedback on protein production at time is often proportional to protein concentration at time , where is known as the delay time. For short delay times, a reaction may be approximated as being instantaneous, and the system as being in quasi equilibrium. However, when the delay is comparable to the characteristic time scales of reaction, the non instantaneous nature of the interactions can no longer be ignored, and delay terms need to be included in the governing equations for the network under study re:mackey95 (); re:bratsun05 ().

Experimental evidence has been mounting that highlights the importance of stochastic effects in transcriptional regulation re:kepler01 (); re:elowitz02 (); re:bratsun05 (); re:maheshri07 (), not only for natural networks, but for engineered gene circuits and networks as well re:elowitz00 (); re:hasty00 (). However, despite the wealth of evidence on the subject, delays in stochastic models of metabolic feedback are very often neglected, possibly because the resulting stochastic equations are no longer Markovian, and hence rarely tractable analytically. Exceptions include the derivation of a two time Fokker-Planck equation and the study of its small delay time limit in re:guillouzic99 (), and results on the bifurcation of the first and second moments of a stochastic linear equation with delay re:mackey95 (); re:lei07 (). We extend these latter results to the analysis of the stationary probability distribution function of a nonlinear model, and discuss in detail the stability of the solutions that results from the interplay of delay and stochasticity.

We focus on a canonical form of a nonlinear Langevin equation with multiplicative or parametric noise and delayed feedback

(1) |

where the constant plays the role of the control parameter, is the intensity of a feedback loop of delay , and is a white, Gaussian noise of intensity . The initial condition is a function specified on . We study the stationary probability distribution function for a range of values of and .

The bifurcation diagram corresponding to Eq. (1) in the absence of noise is known (see, e.g. re:mackey95 ()). Linearization around shows that trajectories decay asymptotically to zero if

(2) |

The boundary separating exponentially decaying solutions from exponentially growing solutions is shown as the solid line in Fig. (4). The upper branch, (I), is defined by and corresponds to a direct bifurcation (real eigenvalue), whereas the lower branch, (II), corresponds to a Hopf bifurcation (complex eigenvalue). In both cases, we show in the figure the case of , and hence the lower brach corresponds to in Eq. (2). The cusp at the intersection of both boundaries is located at .

The stochastic bifurcation analysis of Eq. (1) without delay () is also known re:graham82b (); re:drolet98 (); re:sanmiguel99 (). Analysis of the linearized equation leads to the unphysical conclusion that the bifurcation threshold depends on the order of the statistical moment considered. With the saturating nonlinearity in Eq. (1), stationary probability distributions of can be obtained both below and above threshold, and the location of the threshold properly determined. The stationary distribution function of with is re:graham82b ()

(3) | |||||

(4) |

where the exponent , and is a normalization constant. The solution (4) exists but is not normalizable for and hence it is not a physically admissible solution. Therefore the bifurcation threshold is located at where changes from a delta function at the the origin to a power law at small with an exponential cut off at large . In this case, the existence of parametric fluctuations has no effect on the location of the bifurcation threshold: Both deterministic and stochastic equations exhibit a bifurcation at . In , is unimodal with a divergence at , whereas for the distribution is bimodal reflecting saturation of .

We now turn to the case of delay, . Analytical results for the stability of the trivial solution of the linearization of Eq. (1) have been given in re:mackey95 (); re:lei07 (), and are shown in Fig. (4). The bifurcation threshold of the first moment is shifted relative to the deterministic limit, only bounds have been given for the second moment re:lei07 (), and no results are available for . Given the anomalous behavior described above for the stochastic bifurcation of the linear equation with , it is of interest to determine for the full model of Eq. (1). Unfortunately, the non Markovian character of this equation has precluded progress along these lines re:guillouzic99 ().

We have first extended an existing high order algorithm for the integration of stochastic differential equations re:fox91 () to the case of delay. The algorithm needs to take into account trajectories into the past for an interval , and also new contributions from the stochastic terms that result from the coupling to the delayed feedback. In the numerical results shown below we employ, for technical reasons, an Ornstein-Uhlenbeck stochastic process with intensity and correlation time , the same as the time step in the discretization of Eq. (1). The initial condition considered in all our calculations is a white and Gaussian random process in of zero mean and intensity 1. The time step used is the numerical integration is .

A qualitative view of the bifurcation is given in Fig. 1, which shows a histogram of once trajectories have reached a statistical steady state. For , the histogram is approximately a delta function at . At a critical value , the bifurcation point, a broad distribution emerges, although the most likely value remains . At larger values of , the histogram becomes bimodal. This histogram shown corresponds to the direct bifurcation branch, but a similar graph is obtained for the Hopf bifurcation. Our results for the distribution function in these three ranges of values of are shown in Fig. 2. For is approximately a power law with effective exponent , but with a growing amplitude of at (not shown in the figure). Because of normalization, this growth implies a decaying amplitude for finite , signaling a long transient leading to the build up of the delta function at . Interestingly, the effective power law in the figure , indicating that would not be normalizable. For we do obtain a time independent distribution with . This distribution is normalizable, and represents the stationary distribution above threshold. We finally show in the range of for which the distribution is bimodal. The probability of the most likely value is approximately constant, but we still observe some transients in the vicinity of . Figure 3 shows the results of a power law fit to as a function of the control parameter . We observe a smooth variation of the exponent with that allows a convenient determination of , the value of for which . This is the method that we have used to determine the bifurcation threshold in all the results presented below.

We summarize our results for the bifurcation diagram of Eq. (1)in Fig. 4. The analytic results for the threshold without noise () are shown for reference, as well as the threshold of of the linearized equation re:mackey95 (), and our numerical estimate. Except in the vicinity of the multi critical point our results are in excellent agreement with the analytic calculation of the linear equation, thus validating the accuracy of the numerical algorithm. As one approaches the point the Hopf frequency approaches zero and it is necessary to integrate the differential equation up to very long times to differentiate between an unstable trivial solution or an oscillation with a very long period. Since bounds on the threshold from for the linearized equation have been given in the literature re:lei07 (), we also show in Fig. 4 the threshold for this moment obtained numerically. Our numerical results do agree with the known threshold for the special point of no delay given in re:graham82b (). Interestingly, the thresholds for and converge to the same values when one moves away from the point , both in the direct and Hopf bifurcation branches. The figure also presents our results for bifurcation threshold defined directly from the stationary probability distribution function as discussed above. Our conclusion is that the stochastic threshold is shifted away from the deterministic threshold except in the special point of no delay (), both for the direct and Hopf bifurcations. This threshold also agrees with that of the first and second moments of the linearized equation in the range of parameters in which both agree. Figure 5 shows the dependence of the shift in threshold as a function of the noise intensity . In analogy with the case of no delay, we find a linear dependence in .

We next turn to the equation for that follows from Eq. (1). The difficulty in obtaining a closed, analytic expression lies in the need to find the joint probability distribution . When is larger than the correlation time of , one can assume statistical independence between and , or , an approximation that has been considered in the literature (e.g., ref. re:bratsun05 ()). However, this assumption does not hold near a bifurcation since characteristic correlation times diverge. Instead we write , where is the conditional probability of finding at given at time . We have then derived a the Fokker-Planck equation for that requires only the determination of , a correlation which we have not been able to compute analytically. We can, however, examine its behavior numerically. We have found that it reaches a stationary function (independent of ), and that for small is well described by (for ). With this empirical relation and some straightforward algebra, we recover the same solution for given in Eq. (4) but with . This solution closely agrees with our numerical results of versus as shown by the solid line in Fig. 3, and with the value of the shift in as a function of the noise intensity shown in Fig. 5.

In summary, many of the qualitative features of stochastic bifurcations under multiplicative noise are preserved under the addition of a delayed feedback. First, the bifurcation remains sharp. Since it is commonly assumed in the literature that the correlation time of , the delay term in Eq. (1)) would effectively act as an additive source of noise, leading perhaps to an imperfect bifurcation. We have shown this not to be the case. Second, and in agreement with the case of no delay (), we observe that the moments of the linearized equation bifurcate at different values of the control parameter, which are themselves different from the threshold predicted from the distribution function of the full nonlinear equation. In contrast with the case of , however, the existence of delay introduces a shift in the bifurcation threshold, both for direct and Hopf bifurcations, shift that goes to zero as . The magnitude of the shift scales linearly with the noise intensity . Finally, we have empirically derived the stationary solution of the distribution that agrees with our numerical determination of the bifurcation threshold. In view of our results, care must be exercised when analyzing bifurcation thresholds in numerical simulations of model gene regulatory networks when the analysis is based on the calculation of moments.

This research has been supported by NSERC Canada.

## References

- (1) R. Driver, Ordinary and Delay Differential Equations, Vol. 20 of Applied Mathematical Sciences (Springer, New York, 1977).
- (2) M. Mackey and L. Glass, Science 197, 287 (1977).
- (3) L. Glass and M. Mackey, From Clocks to Chaos: The Rythms of Life (Princeton University Press, Princeton, NJ, 1988).
- (4) J. Hasty, D. McMillen, F. Isaacs, and J. Collins, Nat. Rev. Gen. 2, 268 (2001).
- (5) M. C. Mackey and I. G. Nechaeva, Phys. Rev. E 52, 3366 (1995).
- (6) D. Bratsun, D. Volfson, L. S. S. Tsimring, and J. Hasty, Proc. Natl. Acad. Sci. USA 102, 14593 (2005).
- (7) M. Elowitz, A. Levine, E. Siggia, and P. Swain, Science 297, 1183 (2002).
- (8) N. Maheshri and E. K. O’Shea, Annu. Rev. Biophys. Biomol. Struct. 36, 413 (2007).
- (9) T. Kepler and T. Elston, Biophys. J. 81, 3116 (2001).
- (10) M. B. Elowitz and S. Leibler, Nature 403, 335 (2000).
- (11) J. Hasty, J. Pradines, M. Dolnik, and J. Collins, Proc. Natl. Acad. Sci. USA 97, 2075 (2000).
- (12) S. Guillouzic, I. L’Heureux, and A. Longtin, Phys. Rev. E 59, 3970 (1999).
- (13) J. Lei and M. C. Mackey, SIAM Journal of Applied Mathematics 67, 387 (2007).
- (14) R. Graham and A. Schenzle, Phys. Rev. A 25, 1731 (1982).
- (15) F. Drolet and J. Viñals, Phys. Rev. E 57, 5036 (1998).
- (16) M. San Miguel and R. Toral, in Instabilities and Nonequilibrium Structures, V, edited by E. Tirapegui and W. Zeller (Kluwer Academic, The Netherlands, 1999).
- (17) R. F. Fox, Phys. Rev. A 43, 2649 (1991).