Pitchfork and Hopf bifurcation thresholds in stochastic equations with delayed feedback
The bifurcation diagram of a model stochastic differential equation with delayed feedback is presented. We are motivated by recent research on stochastic effects in models of transcriptional gene regulation. We start from the normal form for a pitchfork bifurcation, and add multiplicative or parametric noise and linear delayed feedback. The latter is sufficient to originate a Hopf bifurcation in that region of parameters in which there is a sufficiently strong negative feedback. We find a sharp bifurcation in parameter space, and define the threshold as the point in which the stationary distribution function changes from a delta function at the trivial state to at small (with exactly at threshold). We find that the bifurcation threshold is shifted by fluctuations relative to the deterministic limit by an amount that scales linearly with the noise intensity. Analytic calculations of the bifurcation threshold are also presented in the limit of small delay that compare quite favorably with the numerical solutions even for .
pacs:02.30.Ks, 05.10.Gg, 05.70.Ln, 87.16.Yc, 87.18.Cf
We study the bifurcation diagram of a nonlinear stochastic differential equation that includes delayed feedback. The model equation considered exhibits both pitchfork and Hopf bifurcations. The bifurcation thresholds are obtained as a function of the model parameters, and our 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) re:belair89 (). 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 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 hence 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 pointing to the importance of stochasticity in feedback regulation, 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 linear stochastic 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 here on the normal form for a pitchfork bifurcation augmented with multiplicative or parametric noise (additive noise has been shown to have no effect in the case of a delayed equation re:mackey95 ()), and linear delayed feedback. The latter is sufficient to originate a Hopf bifurcation in some region of parameters (sufficiently strong negative feedback). Our model for the dynamical variable is
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 generally focus on the stationary probability distribution function for a range of values of and .
The bifurcation diagram for the differential delay equation resulting from Eq. (1) with is known (see, e.g. re:mackey95 ()). Linearization around shows that trajectories decay asymptotically to zero according to
The boundary separating exponentially decaying solutions from exponentially growing solutions is shown as the solid line in Fig. (4). The upper branch corresponds to a pitchfork bifurcation (real eigenvalue, Eq. (2)), whereas the lower branch corresponds to a Hopf bifurcation (complex eigenvalue, Eq. (3)). In both cases, we show in the figure the case of . The cusp at the intersection of both boundaries is located at .
The bifurcation threshold of Eq. (1) without delay () is also known, and has been given in re:graham82b (); re:drolet98 (); re:sanmiguel99 (). Recall that an analysis of the linearized equation leads to the unphysical conclusion that the bifurcation threshold depends on the order of the statistical moment of considered. On the other hand, with a saturating nonlinearity in Eq. (1), a stationary probability distribution of exists both below and above threshold, thus allowing a proper determination of the bifurcation point. The stationary distribution function of with has been found to be re:graham82b ()
where the exponent , and is a normalization constant. The solution (5) exists but is not normalizable for and hence it is not a physically admissible solution. Therefore the stochastic bifurcation threshold is located at , point at which changes from a delta function at the the origin to a power law at small with an exponential cut off at large . Contrary to what an analysis of the moments of would indicate, the existence of parametric fluctuations has no effect on the location of the bifurcation threshold: both deterministic and stochastic equations exhibit a pitchfork bifurcation at . We further note that in , is unimodal with a singularity at , whereas for the distribution is bimodal reflecting nonlinear saturation of .
Ii Stochastic bifurcation
We now turn to the case of delay, . Analytical results for the stability of the trivial solution for 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 delay equation result; only bounds can be given for the stability of 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 second 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. Derivation of the algorithm is presented in appendix A. In the numerical calculations to be presented below, the initial condition is a constant function in for each trajectory, with the constant being drawn from a Gaussian distribution of zero mean and variance . The time step used in the numerical integration is . We also present approximate analytic calculations in the limit of small delay time , following the approach of Frank re:frank05 (). Reasonable agreement is found with our numerical calculations for .
ii.1 Pitchfork bifurcation
A qualitative view of the bifurcation of Eq. (1) is given in Fig. (1). Equation (1) has been integrated numerically, and histograms of computed once trajectories approach a statistical steady state. The histograms shown correspond to independent trajectories for each value of . For , the histogram is approximately a delta function at . As discussed further below, we observe long transients in until it eventually decays to . 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. The histogram shown in the figure corresponds to the direct bifurcation branch, but a qualitatively similar graph is obtained for the Hopf bifurcation.
Our results for the stationary distribution function are shown in Fig. (2). For , we expect . We observe instead a very long lasting transient with approximately characterized by a power law distribution, with an effective exponent (leading to a non normalizable distribution). The amplitude of the point (not shown in the figure) grows with time, signaling the build up of the delta function at the origin. Because of overall normalization, growth at implies a decaying amplitude for finite , as shown in the figure. For , we do obtain a time independent power law distribution with exponent . 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 function around the most likely value is approximately constant over time, but we still observe some transients in the region around . Figure (3) shows our results for the exponent of a power law fit to at small , as a function of the value 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 earlier results for the threshold of of the linearized equation re:mackey95 (). Our numerical results do agree with the known threshold for the special point of no delay given in re:graham82b (). The figure presents our results for the stochastic bifurcation threshold defined directly from the stationary probability distribution function as discussed above. We conclude that the stochastic threshold is shifted relative to the deterministic threshold except in the special point of no delay (). We have also examined in Fig. (5) the dependence of the threshold shift as a function of the noise intensity . In analogy with the case of no delay, we find a linear dependence of with .
ii.2 Hopf bifurcation
The calculation of just shown for the case of a pitchfork bifurcation has been repeated in the vicinity of the deterministic Hopf branch shown in Fig. (4). In the deterministic case, the bifurcation leads to oscillation. When fluctuations are added, oscillation amplitudes fluctuate as well. We also observe in this range of parameters a sharp bifurcation threshold which is shifted relative to the deterministic Hopf bifurcation. The bifurcation threshold is verified with the probability distribution function of the maximum amplitude of the Fourier transform of the trajectories. For each trajectory, we calculate the Fourier transform of the stationary solution over a finite window in time, identify its maximum amplitude, and construct an histogram over those maxima along the trajectory and over the ensemble.
Our results for the stationary probability distribution function of the maximum amplitude of the Fourier transform are shown in Fig. (6) for three increasing values of the control parameter . For we observe an effective power law with exponent and hence not a physical distribution. As was the case above, this is manifested by a transient distribution that decreases with time. As the value of is increased, a power law distribution is found with exponent . The distribution obtained is integrable and stationary. For yet larger values of , the distribution becomes bimodal.
Figure (7) shows the results of a power law fit to the resulting distributions at small for a range of values of of the two probability distribution functions introduced. We have undertaken this analysis for a range of values of , and the resulting Hopf branch of the bifurcation diagram is shown in Fig. (4). Note that it is also shifted relative to the branch corresponding to the deterministic equation.
Iii Fokker-Planck equation
We next turn to an approximate analytic calculation of the stationary distribution function for Eq. (1). The difficulty in obtaining a closed, analytic expression lies in the non Markovian nature of Eq. (1) and the associated need to find the joint probability distribution . When is larger than the correlation time of , the derivation is simplified by assuming statistical independence between and , or . This approximation has already been considered in the literature (e.g., ref. re:bratsun05 ()). However, the assumption of independence does not hold near a bifurcation since the characteristic correlation time diverges re:berbert09 (). We instead proceed as follows: Define the probability distribution of as . By using the properties of the Dirac delta function, and Eq. (1) one finds
where we have introduced the dummy variables and . The second term in the right hand side of Eq. (6) can be written as
by using the Furutsu-Novikov theorem. By using the properties of the delta function, one can write
and since from Eq. (1), we find
where we have also used
Thus the second term in the right had side of Eq. (6) is
Consider now the first term in the right hand side of Eq. (6). By using the law of total expectation,
We anticipate that independence, , does not hold near the bifurcation threshold. We write instead , where is the conditional probability of finding at given the information that . Then
The last integral represents the conditional expected value of given ,
We have not been able to determine analytically. However, it is possible to derive an approximate expression under the assumption that the delay is small re:frank05 (). We write drift terms in Eq. (1) in the Ito representation (recall to we were working with the Stratonovich representation), re:gardiner85 (). Define , , and , so that . Furthermore, since we are in the small time delay approximation, one can assume that the zero-th order approximation of the stationary conditional probability distribution function is Gaussian re:risken89 () and given by,
The approximate drift term of the Fokker-Planck equation is then,
Integrating, one finds,
with stationary solution,
This stationary distribution is found to agree quite well with our numerical determination of (Fig. (2)) around the pitchfork branch, but fails around the Hopf branch. Furthermore, we can now introduce an analytic determination of the bifurcation threshold as the point in which (Eq. (22)) ceases be normalizable: . The bifurcation threshold is located at
This prediction is also in very good agreement with our numerical results for the pitchfork bifurcation branch.
We have also verified our results for the conditional average with a direct numerical integration of the model equation. We show our results for as a function of in Fig. (8). Equation (1) is integrated numerically until it reaches a statistical stationary state. For every value of the dynamical variable , the average of its value at time is collected for independent trajectories in a time window . We note that for large values of , the results become less accurate because of fewer data points in this region. This analysis has been repeated for a range values of , , and . As shown in the figure, we observe good agreement between the numerical results and Eq. (20) in the region around .
In summary, we find that the bifurcation point in our stochastic differential equation with delay remains sharp. This is not a straightforward observation since the delay term in Eq. (1) could effectively act as an additive source of noise, and lead to an imperfect bifurcation instead. This does not appear to be the case. We also note that long transients can be expected below the bifurcation threshold, as all trajectories eventually decay to the trivial solution . Second, we observe that different moments of obtained from the linearized equation bifurcate at different values of the control parameter. When a saturating nonlinearity is introduced into the model, all moments bifurcate at . By defining the bifurcation point in the stochastic case as that in which the power law form of becomes normalizable (), we show that the bifurcation threshold is shifted relative to that of the underlying deterministic equations. The shift goes to zero as , and otherwise it scales linearly with the noise intensity . We have also derived and approximate expression for the stationary distribution function in the limit of small delay time . The threshold location obtained from this approximation agrees well with our numerical determination even when .
This research has been supported by NSERC Canada. MG acknowledges funding by the FQRNT. Computational resources have been provided by CLUMEQ.
- (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. Belair and M. Mackey, J. Dyn. Diff. Eq. 1, 299 (1989), journal of Dynamics and Differential equations.
- (5) J. Hasty, D. McMillen, F. Isaacs, and J. Collins, Nat. Rev. Gen. 2, 268 (2001).
- (6) M. C. Mackey and I. G. Nechaeva, Phys. Rev. E 52, 3366 (1995).
- (7) D. Bratsun, D. Volfson, L. S. S. Tsimring, and J. Hasty, Proc. Natl. Acad. Sci. USA 102, 14593 (2005).
- (8) T. Kepler and T. Elston, Biophys. J. 81, 3116 (2001).
- (9) M. Elowitz, A. Levine, E. Siggia, and P. Swain, Science 297, 1183 (2002).
- (10) N. Maheshri and E. K. O’Shea, Annu. Rev. Biophys. Biomol. Struct. 36, 413 (2007).
- (11) M. B. Elowitz and S. Leibler, Nature 403, 335 (2000).
- (12) J. Hasty, J. Pradines, M. Dolnik, and J. Collins, Proc. Natl. Acad. Sci. USA 97, 2075 (2000).
- (13) S. Guillouzic, I. L’Heureux, and A. Longtin, Phys. Rev. E 59, 3970 (1999).
- (14) J. Lei and M. C. Mackey, SIAM J. Appl. Math. 67, 387 (2007).
- (15) R. Graham and A. Schenzle, Phys. Rev. A 25, 1731 (1982).
- (16) F. Drolet and J. Viñals, Phys. Rev. E 57, 5036 (1998).
- (17) M. San Miguel and R. Toral, in Instabilities and Nonequilibrium Structures, V, edited by E. Tirapegui and W. Zeller (Kluwer Academic, The Netherlands, 1999).
- (18) R. F. Fox, Phys. Rev. A 43, 2649 (1991).
- (19) T. Frank, Phys. Rev. E 72, 011112 (2005).
- (20) J. Berbert, M. Gaudreault, and J. Viñals, , unpublished.
- (21) C. W. Gardiner, Handbook of Stochastic Methods for Physics, Chemistry and the Narutal Sciences, second edition (Springer-Verlag, Berlin, Heidelberg, 1985).
- (22) H. Risken, The Fokker-Planck Equation : Methods of Solution and Applications (Springer, Berlin, 1989).
Appendix A Algorithm for stochastic differential equations with delay
We summarize in this appendix the numerical algorithm developed to integrate Eq. (1). Formal integration yields,
We now approximate and so that to first order in we have,
The nonlinear term also must be expanded around
In order to calculate the integrals containing the random process , we define
If is a Gaussian process of zero mean, and are also Gaussian variables of zero mean, and correlations
The three remaining integrals can be expressed in terms of and as
The Gaussian variables and can be simulated with two Gaussian random variables, and , of zero mean and variance one:
Combining all our results, we write the iteration of our algorithm,