Pitchfork and Hopf bifurcation thresholds in stochastic equations with delayed feedback

Pitchfork and Hopf bifurcation thresholds in stochastic equations with delayed feedback

Mathieu Gaudreault, Françoise Lépine, and Jorge Viñals Department of Physics, McGill University, Montreal, Quebec, Canada, H3A 2T8.
July 4, 2019
Abstract

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

I Introduction

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

(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 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

(2)
(3)

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 ()

(4)
(5)

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

(6)

where we have introduced the dummy variables and . The second term in the right hand side of Eq. (6) can be written as

(7)

with

(8)

by using the Furutsu-Novikov theorem. By using the properties of the delta function, one can write

(9)

and since from Eq. (1), we find

(10)

where we have also used

(11)

Thus the second term in the right had side of Eq. (6) is

(12)

Consider now the first term in the right hand side of Eq. (6). By using the law of total expectation,

(13)

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

(14)

The last integral represents the conditional expected value of given ,

(15)

Substitution of Eqs. (12) and (14) into Eq. (6) leads to the Fokker-Planck equation,

(16)

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,

(17)

The approximate drift term of the Fokker-Planck equation is then,

(18)

Integrating, one finds,

(19)

An approximate expression for the conditional average of given is obtained by comparing Eq. (19) with Eq. (14),

(20)

Further substitution of Eq. (20) into Eq. (16) results in a closed form for the steady state distribution,

(21)

with stationary solution,

(22)

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

(23)

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 .

Acknowledgments

This research has been supported by NSERC Canada. MG acknowledges funding by the FQRNT. Computational resources have been provided by CLUMEQ.

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. 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,

(24)

with

(25)

We now approximate and so that to first order in we have,

(26)

The nonlinear term also must be expanded around

(27)

We substitute Eqs. (26) and (27) into equation Eq. (24) and find,

In order to calculate the integrals containing the random process , we define

(28)
(29)

If is a Gaussian process of zero mean, and are also Gaussian variables of zero mean, and correlations

(30)
(31)
(32)

The three remaining integrals can be expressed in terms of and as

(33)
(34)
(35)

The Gaussian variables and can be simulated with two Gaussian random variables, and , of zero mean and variance one:

(36)
(37)

Combining all our results, we write the iteration of our algorithm,

(38)

FIGURES

Figure 1: Long time histogram of (in grey scale) as a function of the control parameter with , , and . The histograms have been collected in the time interval and further averaged over independent runs. In the absence of noise, the critical value of the control parameter for instability is . We find instead that the bifurcation from a delta function to a power law distribution occurs at instead for this set of parameters. Fluctuations around are observed for due to the finite length of the time series.
Figure 2: Probability distribution function for , and at the times given, and averaged over independent realizations. Values of the control parameter shown are: (a) with , (b) with , and (c) with . The distributions in (a) show a clear transient, whereas those in (b) and (c) are stationary. The solid line shows the power law at small ; the domain covered by the line indicated the range of data that were used to estimate , and is placed above or below the curves for clarity. The dashed line is our approximate determination of in the limit of small .
Figure 3: Results of a power law fit to for small . The fitted value of the exponent is shown as a function of the control parameter . We define the bifurcation threshold for the stochastic problem when , or for this parameter set (, , and ). The solid line follows from our approximate determination of in the limit of small . There are no adjustable parameters.
Figure 4: Numerically determined bifurcation diagram for Eq. (1) with () defined as the point in the plane for which , the exponent of the stationary probability distribution function. For reference, we also show the exact bifurcation diagram for the deterministic equation (solid line), and the exact results for the bifurcation threshold of the first moment from the linearized equation re:mackey95 () (dashed line). The two points labeled by lie on the line , and are known results for the case of no delay for (from left to right) , and for deterministic case of . The dotted line is the approximate threshold [Eq. (23)], and the labels () are the numerically calculated Hopf branch from the probability distribution function of the maximum amplitude of the Fourier transform of the trajectories.
Figure 5: Bifurcation threshold from as a function of noise intensity for and . Time averages used for the determination of are in , and independent realizations have been considered. The line in the figure is the prediction from our approximate determination of the stationary probability distribution function. There are no adjustable parameters.
Figure 6: Probability distribution function of the maximum amplitude of the Fourier transform for , and at the times given, and averaged over independent realizations. Values of the control parameter shown are: (a) with , (b) with , and (c) with , as shown by the solid lines, placed above the curves for clarity. The solid lines extend over the range used for estimation of . The distributions in (a) show a clear transient, whereas those in (b) and (c) are stationary.
Figure 7: Results of a power law fit to () and to the probability of the maximum amplitude of the Fourier transform (), for small or small . The fitted value of the exponent is shown as a function of the control parameter . We define the bifurcation threshold for the stochastic problem when , or (), and () for this parameter set (, , and ).
Figure 8: Ensemble average of given , . We have set and The average is computed over realizations for a time interval . The dashed line is , the approximate analytic result in the limit of small .
Comments 0
Request Comment
You are adding the first comment!
How to quickly get a good reply:
  • Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
  • Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
  • Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
""
The feedback must be of minimum 40 characters and the title a minimum of 5 characters
   
Add comment
Cancel
Loading ...
33047
This is a comment super asjknd jkasnjk adsnkj
Upvote
Downvote
""
The feedback must be of minumum 40 characters
The feedback must be of minumum 40 characters
Submit
Cancel

You are asking your first question!
How to quickly get a good answer:
  • Keep your question short and to the point
  • Check for grammar or spelling errors.
  • Phrase it like a question
Test
Test description