Enhancing noise-induced switching times in systems with distributed delays
EEnhancing noise-induced switching times in systems withdistributed delays
Y.N. Kyrychko ∗ Department of Mathematics, University of Sussex,Falmer, Brighton, BN1 9QH, United Kingdom
I.B. Schwartz
US Naval Research Laboratory, Code 6792,Nonlinear System Dynamics Section,Plasma Physics Division, Washington, DC 20375, USA (Dated: today)
Abstract
The paper addresses the problem of calculating the noise-induced switching rates in systems withdelay-distributed kernels and Gaussian noise. A general variational formulation for the switchingrate is derived for any distribution kernel, and the obtained equations of motion and boundaryconditions represent the most probable, or optimal, path, which maximizes the probability ofescape. Explicit analytical results for the switching rates for small mean time delays are obtainedfor the uniform and bi-modal (or two-peak) distributions. They suggest that increasing the widthof the distribution leads to an increase in the switching times even for longer values of mean timedelays for both examples of the distribution kernel, and the increase is higher in the case of thetwo-peak distribution. Analytical predictions are compared to the direct numerical simulations,and show excellent agreement between theory and numerical experiment. ∗ [email protected] a r X i v : . [ phy s i c s . d a t a - a n ] M a y any real world dynamical systems exhibit complex behavior often inducedby intrinsic time delays. In addition, the majority of such systems are influencedby both external and internal random perturbations. An important practicalproblem, therefore, is to understand how random disturbances are organizedsuch that the dynamics escape from a stable attractor and exhibit new behav-ior. Although the noise amplitudes are small, the resulting change is a largefluctuation out of the basin of attraction. In this paper, we study the influenceof noise-induced large fluctuations on dynamical systems, where the time delayis not taken as a constant but is rather chosen from a given distribution. Weuse a variational approach to calculate the switching rates out of the basin ofattraction of the stable equilibrium for a general kernel of the delay distribution,as well as general nonlinearity. Taking two particular commonly-used examplesof the distribution kernel, namely, the uniform and bi-modal kernels, we ana-lyze how the width of the distribution affects the switching rates. Our resultssuggest that the switching is affected not only by the mean time delay, but alsoby the width of the delay distribution. Specifically, if the distribution width isincreased, this leads to an increase in the switching times. I. INTRODUCTION
Time-delayed models have been extensively applied in various disciplines, including neu-roscience, quantum optics, laser dynamics, mechanical and chemical oscillators, populationdynamics, mathematical epidemiology, and many others [1]. In these systems, time delaysare often used to account for a finite propagation time, the time it takes for a signal to beprocessed and fed back, and significant time delays can arise due to the separation betweenthe subsystems. In neural systems, delays represent the time it takes for neural signalsto propagate and be processed [2, 3], and in semiconductor lasers the delay is determinedby the propagation path, and is much longer than the internal oscillation periods [4]. Inquantum networks, in order to model non-Markovian aspects of quantum dynamics, one hasto account for delayed interactions between network nodes when exchanging photons, andtime-delayed feedback control can be used to create and stabilize entangled states [5, 6].In order to describe population fluctuations in ecological setting, it is important to include2aturation time, the time it takes to develop into a reproducing adult [7–9]. In epidemiol-ogy, time delays can represent latency or temporary immunity periods affecting the spreadof the infection [9–11], and in gene regulatory networks they arise during the transcriptionand translation of mRNA [12–14]. Finally, robotic systems such as swarms and their controlhave attracting states that only can be captured with delayed communication and controlactuation [15, 16].In addition to intrinsic time delays, majority of real systems are also affected by bothexternal and internal random perturbations, which necessitate the inclusion of noise intothe dynamical models, and hence, stochastic differential delay equations (DDEs) are used toanalyze a wide range of systems. For instance, gene regulatory networks are often modeledas stochastic birth-death processes with time delays [12], and anomalies during El Ni˜nohave been analyzed using delayed-oscillator models, excited by external weather noise [17–19]. Recently, time-delayed feedback control has been applied to noise-induced chimerastates in FitzHugh-Nagumo networks with non-local coupling, and it has been shown thatfor some values of time delays, it can lead to the so-called period-two coherence resonancechimera [20]. Much of the analysis of stochastic DDEs involving local fluctuations has beendone using delayed Fokker-Plank equations, and several small delay approximation methodshave been suggested. Furthermore, a perturbation theory method has been developed todetermine single time point probability densities and mean delays, as well as autocorrelationfunctions in stochastic SDEs in the context of delayed Fokker-Plank equations [21, 22].In this paper we focus on understanding large fluctuations in dynamical systems withdistributed time delays, which have one or more stable steady states. In the absence of noise,the solutions will go to a stable steady state; however, when the noise is included, even ifsmall in intensity, it leads to fluctuations around the stable steady state. Furthermore, noisecan force the system out of the stable steady state, or it can lead to switching to anotherstable steady state. This phenomenon is similar to the tipping process in dynamical systems,where relatively small changes in input can lead to sudden and disproportional changes inoutput, for example, due to a slow variation in parameter (B-tipping) or inclusion of noise(N-tipping) [23–25]. For systems not in thermal equilibrium, as in the case of noise-inducedswitching between states, the probability distribution is no longer of the Boltzmann form.Many results have been obtained for dynamical systems without delay driven by whiteGaussian noise and for Markovian reaction and population systems, cf. [26–37].3n the other hand, noise-induced switching has been studied for bistable systems inthe presence of time delay to demonstrate residence times depending on the noise level, aswell as to analytically calculate the expressions for the autocorrelation function and powerspectrum [38, 39]. Switching entire behaviors for complex discrete delay-coupled swarmingsystems has been observed in [40]. Noise-induced switching rates between stable steadystates, and the rates of noise-induced extinction in the case of constant time delays, haverecently been analyzed in general in [41]. In this paper we will study a time-delayed systemdriven by noise, where the time delay is given by an integral with a memory kernel in theform of the prescribed delay distribution function. The research in this paper is motivatedby the effects of distributed delays on the system dynamics in various settings; for example,in models coupled oscillators [42–44], in traffic dynamics to describe the memory effects ondrivers [45], to represent long-range interactions between neurons [3, 46], modeling waitingtimes in epidemiological models [47], as well as maturation periods in population dynamicsmodeling [48, 49].The outline of the paper is as follows. In Section II we apply variational approach toformulate the problem in order to calculate the rate of switching out of the basin of attractionof the steady state for a general delay distribution and general nonlinearity. In Section III,we consider the case of the uniformly distributed delay kernel, and show how the width ofthe distribution influences the switching rate. The case of the two-peak delay distribution isanalyzed in Section IV, and the results are compared to the case of the uniform distribution.For both types of distributions, we perform numerical simulations, which are presented inSection V, and the agreement between theory and analytical results is excellent.
II. GENERAL CASE OF THE DELAY DISTRIBUTION KERNEL
We consider a switching process in a one-dimensional system with distributed time delayand noise, which has the form˙ x ( t ) = f (cid:18) x ( t ) , (cid:90) ∞ g ( σ ) x ( t − σ ) dσ (cid:19) + η ( t ) , (1)where the distribution kernel g ( · ) is assumed to be non-negative and normalized to unity,i.e. g ( s ) ≥ , (cid:90) ∞ g ( s ) ds = 1 , ( t ) is the noise, and f ( · ) is a nonlinear function. If g ( u ) = δ ( u ), one obtains a systemwithout time delays, and if g ( u ) = δ ( u − τ ), this leads to the case of discrete time delay.In the case of Gaussian noise, if D is the noise intensity, one can characterize noiserealizations η ( t ) by its probability density functional P η [ η ( t )] ≈ exp( −R η /D ) with R η [ η ( t )] = 14 (cid:90) ∞−∞ dtdt (cid:48) η ( t ) ˆ F ( t − t (cid:48) ) η ( t (cid:48) ) , (2)where ˆ F ( t − t (cid:48) ) / D is the inverse of the pair correlator of η ( t ) [50]. Noise intensity D isassumed to be sufficiently small, so that sample paths will limit on an optimal path as D → x A , and a saddle point, x S , satisfying f ( x A , x A ) = f ( x S , x S ) ≡
0. For thetopological switching setup, we suppose the saddle lies on the boundary of the basin ofattraction of x A . That is, we will consider starting initially in the basin of attraction ofattractor x A , and consider the dynamics of switching as escape from the basin of attraction,which we define as a large fluctuation. In the asymptotic limit as the noise intensity goes tozero, the path away from x A will pass through the saddle point, x S .When compared to the effective barrier height, we assume the noise is small, and the meantime to switch is much longer than the relaxation time. Thus, the occurrence of switchingis expected to be a rare event. Assuming the event lies in the tail of the distribution, theprobability of a switching is an exponential distribution, given as P x [ x ] ∝ exp( − R/D ) , R = min R [ x, η, λ ] , (3)where R [ x, η, λ ] = R η [ η ( t )] + (cid:90) ∞−∞ λ ( t ) (cid:20) ˙ x ( t ) − f (cid:18) x, (cid:90) ∞ g ( σ ) x ( t − σ ) dσ (cid:19) − η ( t ) (cid:21) dt, (4)where λ ( t ) is a time-dependent Lagrange multiplier. Note that the noise source η ( t ) may befrom any continuous or discrete distribution. Here we assume uncorrelated Gaussian noisedefined in equation (2) so that we substitute it into equation (4) as R η [ η ( t )] = 14 (cid:90) η ( t ) dt. (5)The main goal to determine the probability of switching is to compute the exponent R ,which, for reasons that will be made clear below, we define as the action, similar to theaction in classical mechanics. 5 . The Variational Equations of Motion To simplify the analysis, we define the following convolution:( g (cid:126) x )( t ) ≡ (cid:90) ∞ g ( σ ) x ( t − σ ) dσ. (6)In order to find the exponent R in equations (4) and (5), we look for equations that describethe maximum probability of reaching the state x S from the initial state x A . The variation( δ R ) is obtained by varying deviations from the path that minimize R . Variation withrespect to noise η ( t ) is η ( t ) = 2 λ ( t ).The variation with respect to the Lagrange multiplier λ ( t ) has the form˙ x ( t ) = f ( x ( t ) , ( g (cid:126) x )( t )) + 2 λ ( t ) . Finally, we look at the variation with respect to x , which is given by δ R δx = R [ x + µ, η, λ ] − R [ x, η, λ ]= (cid:90) ∞−∞ λ ( t ) [ ˙ µ ( t ) − f ( x ( t ) + µ ( t ) , ( g (cid:126) x )( t ) + ( g (cid:126) µ )( t ))+ f ( x ( t ) , ( g (cid:126) x )( t ))] dt. Taylor expanding the last expression gives δ R δx = (cid:90) ∞−∞ λ ( t ) (cid:90) ∞ g ( σ ) [ ˙ µ ( t ) − µ ( t ) f ( x ( t ) , ( g (cid:126) x )( t )) (7) − f ( x ( t ) , ( g (cid:126) x )( t )) µ ( t − σ )] dσdt + O ( µ ) , where f i ( · , · ), i = 1 , f = ∂f ( x ( t ) , ( g (cid:126) x )( t )) ∂x and f = ∂f ( x ( t ) , ( g (cid:126) x )( t )) ∂ ( g (cid:126) x )( t ) . Evaluating the second integral in (7) gives the equation of motion for λ in the form: − ˙ λ ( t ) = λ ( t ) f ( x ( t ) , ( g (cid:126) x )( t )) (8)+ (cid:90) ∞ λ ( t + σ ) g ( σ ) f [ x ( t + σ ) , ( g (cid:126) x )( t + σ )] dσ. Note that unlike the equations for the noise and the state functions, the last equationcontains both delayed and advanced terms, making it an acausal equation of motion for λ ( t ). Combining the derived variational results leads to a system of equations in the form6 x ( t ) = f ( x, ( g (cid:126) x )( t )) + 2 λ ( t ) , (9)˙ λ ( t ) = − λ ( t ) f ( x ( t ) , ( g (cid:126) x )( t )) (10) − (cid:90) ∞ λ ( t + σ ) g ( σ ) f [ x ( t + σ ) , ( g (cid:126) x )( t + σ )] dσ, which has the Hamiltonian H ( x, ( g (cid:126) x )( t ) , λ ) = λ + λf ( x ( t ) , ( g (cid:126) x )( t )) . (11) B. Boundary conditions
To solve for the most likely, or optimal, path along which the system switches from theattractor, we assume we start at or near x A , and compute the path to the saddle point x S using the equation (10). In order to compute the actual path, boundary conditions arerequired on the real line. We assume that at steady state, the noise goes to zero. Therefore,we have that as t → ∞ , ( x ( t ) , λ ( t )) → ( x S , t → −∞ , we have( x ( t ) , λ ( t )) → ( x A , η ( t ) ≡
0, we can linearize the equation (1) about a steady state, ¯ x :˙ X ( t ) = f (¯ x, ¯ x ) X ( t ) + f (¯ x, ¯ x )( g (cid:126) X )( t ) . (12)In particular, when ¯ x = x A , the spectrum of the characteristic equation has values in theleft half of the complex plane. In contrast, when ¯ x = x S , at least one eigenvalue has positivereal part.On the other hand, since the equation (8) is linear in λ , at the steady state ¯ x we have˙Λ( t ) = − Λ( t ) f (¯ x, ¯ x ) − f (¯ x, ¯ x ) (cid:90) ∞ g ( σ )Λ( t + σ ) dσ. (13)Assuming Λ( t ) = e αt , the characteristic equation for Λ is, − α − f (¯ x, ¯ x ) − f (¯ x, ¯ x )[ {L g } ( − α )] = 0 , (14)where α is an eigenvalue of a Jacobian, and {L g } ( z ) = (cid:90) ∞ e − zu g ( u ) du,
7s the Laplace transform of the function g ( u ). Moreover, it is easy to show that the spectrumof Λ have opposite signs to that of X . Therefore, the attractor becomes a saddle, and saddleremains a saddle. This result will hold true if the mean value of the delay distribution issmall.The interpretation of the local saddle structure and the escape path is such that atattractor x A , the conjugate variable supplies an unstable direction for escape through x S .Escape occurs along the most probable path which connects saddles ( x A ,
0) to ( x S ,
0) asa heteroclinic orbit. That is, as t → −∞ , the path approaches ( x A ,
0) along its unstablemanifold, while as t → ∞ , the path asymptotes to ( x S ,
0) along the stable manifold.
C. Perturbation Theory
Taking into account the above-mentioned assumptions on the equations of motion andasymptotic boundary conditions, the problem of finding the action R is formulated as calcu-lating the solutions to a nonlinear two-point boundary problem. We assume that this solutionexists in the non-delayed case, and when the delay is non-zero, the corresponding solutionremains close to the non-delayed solution. Thus the variational problem of finding the ac-tion is formulated as finding its perturbation to non-delayed case. We assume that if x ( t )is a solution for a non-delayed problem, then the perturbations (cid:20) x ( t ) − (cid:90) ∞ g ( z ) x ( t − z ) dz (cid:21) should remain small. The action can be written as the following perturbation problem R [ x, η, λ ] = R [ x, η, λ ] + R [ x, η, λ ] , where R [ x, η, λ ] = 14 (cid:90) η ( t ) dt + (cid:90) λ ( t )[ ˙ x ( t ) − f ( x, x ) − η ( t )] dt, (15)and R [ x, η, λ ] = (cid:90) λ ( t ) [ f ( x ( t ) , x ( t )) − f ( x ( t ) , ( g (cid:126) x )( t ))] dt. (16)The minimizing solution can be found by first analyzing the equations that minimize R ,which are denoted by [ x o , η o , λ o ], and then evaluating the first order correction at the zerothorder solution; i.e., R [ x , η , λ ].From equations (15) and (11), it is easy to see that the the optimal path for the non-delayed case is given by λ = − f ( x , x ), implying that ˙ x = − f ( x , x ), which is just8 time-reversed trajectory from the zero delay case. The optimal noise is then given by η ( t ) = 2 λ ( t ).In order to make further analytical progress, we need to specify a particular choice ofthe distribution kernel g ( u ) and of the function f ( · ). In this paper we consider two partic-ular distributions, namely, a uniform and a two-peak distribution kernel, and a quadraticnonlinearity in the function f ( · ). III. UNIFORMLY DISTRIBUTED DELAY
In this section we consider the case of the uniformly distributed kernel that can be writtenas follows, g ( u ) = ρ , τ − ρ ≤ u ≤ τ + ρ, , otherwise. (17)With this distribution kernel, we choose the function f ( · ) to be f ( x ( t ) , ( g (cid:126) x )( t )) = x (1 − x ) − γ ρ (cid:90) τ + ρτ − ρ x ( t − z ) dz. (18)Substituting the expression (18) into the right-hand side of the equation (1), one obtains˙ x ( t ) = x (1 − x ) − γ ρ (cid:90) τ + ρτ − ρ x ( t − z ) dz + η ( t ) . (19) A. Steady State Stability
In the absence of noise, the equation (19) has two steady states, x A = 1 − γ and x S = 0 , ≤ γ ≤ . To study the stability of these steady states, we substitute the Laplace transform of thedistribution kernel {L g } ( − α ) = 12 ρα e ατ ( e αρ − e − αρ ) = e ατ sinh( αρ ) αρ . into the characteristic equation (14) for the conjugate variable, Λ. (Note that since thespectra for X and Λ characteristic equations are of opposite sign, we only examine the one9or Λ.) Using the Taylor expansion for small mean time delays (and hence, for small ρ ), wecan simplify the characteristic equation and find eigenvalues α as − α − x + γ (1 + ατ ) = 0 = ⇒ α = 1 − x − γγτ − . At the steady state ¯ x = x A , we have α = − γ − γτ > α = − γ − γτ < X .Similar but opposite sign results hold for x S .Therefore, when t → −∞ , the solution of the linearized problem tends to ( x A , x S ,
0) as t → ∞ . B. Perturbation of the Action
In order to find the optimal path from x A to x S , which minimizes R , we use the factthat the optimal path solutions satisfy ˙ x o ( t ) = − f ( x o ( t ) , x o ( t )) and λ o = − f ( x o ( t ) , x o ( t )).For the vector field without delay, we find that x o = x A e x A t , (20)that as t → ∞ , x o ( t ) →
0, and as t → −∞ , x o ( t ) → x A and R ( τ, ρ ) = x A . (21)Since f ( x ( t ) , x ( t )) − f ( x ( t ) , ( g (cid:126) x )( t )) = − γ (cid:20) x ( t ) − ρ (cid:90) τ + ρτ − ρ x ( t − z ) dz (cid:21) , the first order action can be found as R ( τ, ρ ) = − γ (cid:90) ∞−∞ ˙ x o ( t ) (cid:20) x o ( t ) − ρ (cid:90) τ + ρτ − ρ x o ( t − z ) dz (cid:21) dt. (22)= − γ (cid:90) ∞−∞ − x o ( x A − x o ) (cid:20) x o ( t ) − x A − ρ ln (cid:18) e x A ( t − τ − ρ ) e x A ( t − τ + ρ ) (cid:19)(cid:21) dt. Combining the above calculations, the expression for R ( τ, ρ ) becomes R ( τ, ρ ) = x A − γx A − γρ · x A ( τ + ρ ) e x A ( τ + ρ ) − γρ · x A ( τ − ρ ) e x A ( τ − ρ ) − , (23)where we have an explicit dependence on the mean time delay τ and the distribution width10 . For small distribution widths ρ , this expression can be further expanded as R ( τ, ρ ) ≈ x A − γ (cid:18) − x A x A e x A τ ( e x A τ − x A τ − e x A τ − (cid:19) (24)+ 13 x A e x A τ ( − e x A τ + 3 + x A τ e x A τ + 4 x A τ e x A τ + x A τ )( e x A τ − ρ + O ( ρ ) . If ρ = 0, then we recover the same expression as in [41], which was derived for the case ofthe discrete time delay. IV. A TWO-PEAK DELAY DISTRIBUTION KERNEL
In this section we analyze the effect of distributed delay on the switching rates by consid-ering a distribution kernel g ( u ) in the form of two discrete time delays with an average delay τ , which are separated by a time interval 2 ρ [51]. Under this approximation, the distributionkernel g ( u ) can be written as g ( u ) = δ ( u − τ − ρ ) + δ ( u − τ + ρ )2 , where δ is the Dirac delta function. For simplicity, we introduce the notation as follows τ + ρ = τ , τ − ρ = τ , and x τ = x ( t − τ ) , x τ = x ( t − τ ) . The calculations for the variation with respect to noise and λ ( t ) will be similar to those inSection II, so we just state the one variation with respect to x ( t ), since the problem nowinvolves multiple delays. Looking at the deviation with respect to x ( t ) gives: δ R δx = − (cid:90) µ ( t ) (cid:20) ˙ λ ( t ) + λ ( t ) δf ( x, x ( t − τ ) , x ( t − τ )) δx ( t ) − λ ( t + τ ) δf ( x ( t + τ ) , x ( t ) , x ( t + 2 ρ )) δx ( t − τ ) − λ ( t + τ ) δf ( x ( t + τ ) , x ( t − ρ ) , x ( t )) δx ( t − τ ) (cid:21) dt. x ( t ) = f ( x ( t ) , x ( t − τ ) , x ( t − τ )) + 2 λ ( t ) , ˙ λ ( t ) = − λ ( t ) δf ( x, x ( t − τ ) , x ( t − τ )) δx ( t ) − λ ( t + τ ) δf ( x ( t + τ ) , x ( t ) , x ( t + 2 ρ )) δx ( t − τ ) − λ ( t + τ ) δf ( x ( t + τ ) , x ( t − ρ ) , x ( t )) δx ( t − τ ) . We remark that since the problem is one with multiple delays, the equations of motion canbe derived using the Hamiltonian: H ( x, x τ , x τ , λ ) = λ + λf ( x, x τ , x τ ) , with the general acausal equations of motion:˙ x o = ∂H∂λ ( x o , x o τ , x o τ , λ o ) , ˙ λ o = − ∂H∂x ( x o , x o τ , x o τ , λ o ) − ∂H∂x τ ( x o ( t + τ ) , x o ( t ) , x o ( t + 2 ρ ) , λ o ( t + τ )) − ∂H∂x τ ( x o ( t + τ ) , x o ( t − ρ ) , x o ( t ) , λ o ( t + τ )) . In order to compute the optimal path, we assume that the solution exists for τ , = 0,and the solutions for τ , (cid:54) = 0 remain close to the non-delayed solution, and as before weassume that perturbations to the optimal path when τ , (cid:54) = 0 remain small. This meansthat if x ( t ) is a solution for τ , = 0, then the perturbations δx τ i ≡ x ( t ) − x ( t − τ i ), i = 1 , R [ x, η, λ ] and R [ x, η, λ ], where R [ x, η, λ ] = 14 (cid:90) η ( t ) dt + (cid:90) λ ( t )[ ˙ x ( t ) − f ( x, x, x ) − η ( t )] dt, and R [ x, η, λ ] = (cid:90) λ ( t )[ f ( x ( t ) , x ( t ) , x ( t )) − f ( x ( t ) , x ( t − τ ) , x ( t − τ ))] dt, so the action now explicitly depends on two constant time delays τ and τ . In order to com-pare the results of the two-peak distribution kernel with the results on uniform distributionobtained in Section III, as well as to the case of the single discrete time delay considered in[41], we consider function f as follows, f ( x, x τ , x τ ) = x (1 − x ) − γ x τ + x τ ) . η ( t ) = 0, there are two steady states x A = 1 − γ and x S = 0 , ≤ γ ≤ . Linearizing near the steady states, using the characteristic equation (14), yields − α − x + γ ατ + 1 + ατ ) = 0 = ⇒ α = 1 − x − γ ( γ/ τ + τ ) − − x − γγτ − , and the stability characteristics of the attractor and saddle are exactly the same as in theprevious section.In the absence of delay, the computation of R is the same as before, and we only needto compute R . Since f ( x, x, x ) − f ( x, x τ , x τ ) = − γ x − x τ ) + ( x − x τ )] , we can find the first order action as R ( τ , τ ) = − γ (cid:90) ∞−∞ ˙ x o ( t )[ x o ( t ) − x o ( t − τ )] dt − γ (cid:90) ∞−∞ ˙ x o ( t )[ x o ( t ) − x o ( t − τ )] dt. = 2 γ (cid:90) ∞−∞ x A (1 + e x A t ) dt − γx A (cid:90) ∞−∞ e x A t (1 + e x A t ) (1 + e x A ( t − τ ) ) dt − γx A (cid:90) ∞−∞ e x A t (1 + e x A t ) (1 + e x A ( t − τ ) ) dt. Evaluating the integrals in the last expression gives the first order action R ( τ , τ ) and thefull action to the first order in δx ( t ) can be found as R ( τ , τ ) ≈ x A − γx A (cid:18) − e x A τ [ − x A τ − e x A τ ]( e x A τ − + e x A τ [ − x A τ − e x A τ ]( e x A τ − (cid:19) , (25)where τ = τ + ρ and τ = τ − ρ . If the width ρ = 0, the time delays τ and τ coincide, andthe expression (25) becomes the same as the one derived in [41] for the case of the singletime delay. V. NUMERICAL SIMULATIONS
In this section we perform numerical simulations of the system (1) with uniform andtwo-peak distribution kernels, given by (19) in order to compare theoretical predictions forswitching times. Following the same methodology as in [41], we calculate the mean switchingtimes T st as the noise η ( t ) takes the system out of the basin of attraction of the steady state13 IG. 1. (a) Sample switching time series in (19) and (b) a probability density of switching obtainedby numerical simulations for the uniformly distributed kernel. Black squares (circles) denote theattractor (saddle point). Parameters are: γ = 0 . √ D = 0 . τ = 0 . ρ = 0 . x A to x < x S . Numerical simulations are run for some time after the trajectory passes x S toensure that the probability of its return to the basin of attraction of the stable steady state x A is exponentially small. Numerical simulations are performed using the Heun or stochasticpredictor-corrector method [52, 53], in which the solution to equation (1) is approximatedusing the following iterative scheme x n +1 = x n + 12 ( K n + K n +1 ) h + √ Dhu, with a step size h , u is N (0 , K n = x n (1 − x n ) − Y n , K n +1 = ¯ x n +1 (1 − ¯ x n +1 ) − Y n +1 , where the approximation to ¯ x n +1 is given by¯ x n +1 ≈ x n + K n h + √ Dhu.
The terms Y n and Y n +1 representing the delay distribution were computed using a compositetrapezoidal rule in the case of uniform distribution, and for the two-peak distribution thesewere given by Y n = 12 ( x n − k + x n − k ) , Y n +1 = 12 ( x n +1 − k + x n +1 − k ) , where k h = τ − ρ and k h = τ + ρ . 14 IG. 2. A comparison of the first order perturbation theory in time delay τ given by (19) and thelog of the mean switching times obtained by numerical simulations for the uniformly distributedkernel for several values of the distribution width ρ . Solid lines are theoretical predictions, andsymbols are numerical simulations. Other parameters are: γ = 0 . √ D = 0 . τ = 0 . Examples of switching time series for several runs are shown in Figs. 1(a) and 3(a) forthe uniform and bi-modal delay distributions, respectively. One can see that the dynamicsfluctuates around a mean value of the attractor for long period of time until at some pointthere is a large change in which the dynamics passes through the saddle at the origin. Theaccompanying probability densities in panel (b) illustrate the fact that most of the time thedynamics are in the neighborhood of the attractor. Although the probability density plotsshow many paths to going from x A to the saddle at the origin, they also contain the mostprobable path to switch.The switching rate is proportional to the probability of large fluctuations P x [ x ] introducedin Section II, and the switching time is the inverse of the switching rate. Using the equation(3), the switching time can be found as T sw = c sw exp( R sw /D ) , R sw is given by the approximation (24) for the uniformly distributed delay kernel,and the expression (25) in the case of the two delay kernel approximation. Since the noise isGaussian, the constant c sw can be found using the Kramer’s theory [54] at zero delay. Thisconstant for non-zero delay corresponds to the approximate vertical shift of the theoreticalresults when T sw is plotted on the logarithmic scale, and is assumed to have weak dependenceon noise and delay. FIG. 3. a) Sample switching time series in two-peak delay distribution and (b) a probability densityof switching obtained by numerical simulations. Black squares (circles) denote the attractor (saddlepoint). Parameters are: γ = 0 . √ D = 0 . τ = 0 . ρ = 0 . In Figs. 2 and 4, the lines represent theoretical approximations as the first order in τ perturbation theory given by the expressions (24) and (25), and circles are the meanvalues of the numerical simulations taken over 2000 simulations. In Fig 2, we show thecomparison between the theoretical results and numerical simulations of the switching timeas the function of mean time delay τ for the case of uniformly distributed kernel for severalvalues of the distribution width ρ . It can be seen that increasing the mean time delay τ forany distribution width ρ leads to the decrease in the switching times T sw , and theoreticalpredictions align well with numerical simulations all the way up to τ = 1. It is also clearthat increasing the distribution width leads to the increase in the switching times even forsmall values of the mean time delay τ , and this increase in switching times in maintainedfor a whole range of admissible values of τ ( τ < ρ ).A similar situation is observed in Fig 4, which shows the comparison between theory andnumerical simulations for the case of the two-peak distributed kernel, given by two discrete16 IG. 4. A comparison of the first order perturbation theory in time delay τ and the log ofthe mean switching times obtained by numerical simulations for two-peak distribution kernel forseveral values of the distribution width ρ . Solid lines are theoretical predictions, and symbols arenumerical simulations. Other parameters are: γ = 0 . √ D = 0 . τ = 0 . time delays separated by a time interval of 2 ρ . Again, analytical results and numericalsimulations are in an excellent agreement for mean time delays up to τ = 1. The fastestswitching times are achieved when the width ρ is equal to zero, corresponding to the case of asingle discrete time delay in model (1), and when the mean time delay τ is large enough. Asthe distribution width ρ becomes larger, the switching time for the trajectory to be pushedout of the basin of attraction of the steady state increases over the whole range of timedelays τ . This clearly shows that in delayed systems, not only the value of the mean timedelay plays a role in determining the escape rates/times, but the width of the distributionalso has an important influence on the overall dynamics, and the escape rates in particular.Moreover, comparing Figs. 2 and 4, it is worth noting, that while in both cases the switchingtimes grow for larger distribution widths, the growth is more pronounced in the case of thetwo-peak distribution. 17 I. DISCUSSION
In this paper we have considered the problem of finding escape rates from the basin ofattraction of the stable steady state in a dynamical system with distributed time delaysand Gaussian noise. In the absence of noise, the system has stable and unstable steadystates, and the solutions initially in the basin of attraction of the stable steady state remainin the basin. When the noise is included, most of the time the solution trajectories willfluctuate around the stable steady state. However, there exist rare instances where noiseacts as a coherent force, resulting in a trajectory that takes the system out of the attractor’sbasin of attraction. Here we have aimed to understand how particular delay distributionsinfluence the escape rates/times as compared to a single discrete time delay situation. Inorder to calculate the escape rate, the problem has been formulated for a general distributionkernel as a variational problem along the optimal path, which maximizes the probability ofescape. Since noise-induced switching occurrences are rare events, the optimal solution tothe variational problem is valid in the tail of the distribution, which is exponential. Wealso note that even small changes in the action lead to exponential changes in switchingprobability and switching times.In order to make further analytical progress, having obtained the variational equations forthe extreme trajectories for a general distribution, we have then considered two particularcases of the distribution kernel, namely, uniform and two time delay (two-peak) distribu-tions. For both of those distributions, we have been able to obtain closed form analyticalexpressions for the escape rates, which explicitly depend on the mean delay and distributionwidth. The exemplified results hold near a trans-critical bifurcation point. It is known thatthe exponent of the probability of escape from the attractor scales linearly with distancefrom the bifurcation point when the delay is zero. In contrast, we note that the scaling ex-ponent of a saddle-node bifurcation scales as χ / , where χ is the distance to the bifurcationpoint. [41]We have performed numerical simulations for both distribution kernels and compared theresults of the theoretical predictions for the escape time to the numerically calculates ones.We have found very good agreement between theoretical prediction based on the small-delayapproximation, and numerical simulations for mean time delay values up to τ = 1. We havealso compared theoretical and numerical results for larger values of time delay τ , but found18hat they start to diverge for larger time delays, suggesting that the linear approximationdoes not work well for large time delays.We have found that the fastest switching times for all values of the time delay are in thecase when the distribution width is equal to zero for both examples of the distribution kernelconsidered. As the distribution width is increased, the switching times are decreasing evenfor the large values of the mean time delay. Moreover, comparing the change in switchingtimes for uniform and two-peak distribution, we have found that while both do increaseescape times, the latter has a more profound effect. This strongly indicates that not onlythe mean time delay plays a role in the switching dynamics, but the width and the choice ofthe distribution are another two important factors that influence exponentially the switchingrates in time-delayed systems.In this paper we have concentrated on the influence of the delay distribution on theswitching rates in the presence of Gaussian noise and have used two particular distributionkernels. The next step would be to analyze the influence of the Gamma distributed kernel(for weak and strong Gamma kernels) and compare the results to the cases of a single discretetime delay studied in [41], and of the uniform and two-peak distributions considered in thispaper. The model studied in this paper has a saddle point and only a single stable steadystate. Another important future direction of this research would be to analyze the influenceof the delay distribution on the switching times in bistable systems. VII. ACKNOWLEDGMENTS
The authors gratefully acknowledge support from the Office of Naval Research undercontract numbers N0001418WX01225 and N0001412WX20083, and support of the NRLBase Research Program N0001412WX30002. [1] T. Erneux, Applied delay differential equations (Springer, Springer, New York, 2009).[2] J. B´elair, S. A. Campbell, and P. van den Driessche, SIAM J. Appl. Math. , 245 (1996).[3] R. Rahman, K. Blyuss, and Y. Kyrychko, SIAM J. Appl. Dyn. Sys. , 2069 (2015).[4] T. Heil, I. Fischer, W. Els¨asser, J. Mulet, and C. Mirasso, Phys. Rev. Lett. , 795 (2001).[5] H. Pichler and P. Zoller, Phys. Rev. Lett. , 093601 (2016).
6] S. M. Hein, F. Schulze, A. Carmele, and A. Knorr, Phys. Rev. A , 052321 (2015).[7] S.-L. Wu, P. Weng, and S. Ruan, Europ. J. Appl. Math. , 61 (2015).[8] S. Gourley, J.-H. So, and J. Wu, J. of Math. Sci. , 5119 (2004).[9] Y. Kyrychko, S. A. Gourley, and M. V. Bartuccelli, SIAM J. Math. Anal. , 1688 (2006).[10] X. Zhong, S. Guo, and M. Peng, Stochastic Anal. Appl. , 1 (2017).[11] L. Bauer, J. Bassett, P. H¨ovel, Y. N. Kyrychko, and K. B. Blyuss, Chaos: An InterdisciplinaryJournal of Nonlinear Science , 114317 (2017).[12] C. Gupta, J. M. L´opez, R. Azencott, M. R. Bennett, K. Josi´c, and W. Ott, J. Chem. Phys. , 204108 (2014).[13] K. Parmar, K. Blyuss, Y. Kyrychko, and S. Hogan, Comp. and Math. Methods in Medicine , 347273 (2015).[14] J. Feng, S. A. Sevier, B. Huang, D. Jia, and H. Levine, Phys. Rev. E , 032408 (2016).[15] K. Szwaykowska, I. B. Schwartz, L. Mier-y Teran Romero, C. R. Heckman, D. Mox, andM. A. Hsieh, Phys. Rev. E , 032307 (2016).[16] J. Hindes, K. Szwaykowska, and I. B. Schwartz, Phys. Rev. E , 032306 (2016).[17] G. Burgers, Climate Dynamics , 521 (1999).[18] C. Ping, L. Ji, H. Li, and M. Fl¨ugel, Physica D , 301 (1996).[19] A. Keane, B. Krauskopf, and C. Postlethwaite, SIAM J. Appl. Dyn. Sys. , 1229 (2015).[20] A. Zakharova, N. Semenova, V. Anishchenko, and E. Sch¨oll, Chaos , 114320 (2017).[21] T. Frank, Phys. Lett. A , 275 (2006).[22] T. Frank, Phys. Lett. A , 1341 (2016).[23] P. Ashwin, S. Wieczorek, R. Vitolo, and P. Cox, Phil. Trans. R. Soc. A , 1166 (2012).[24] T. Nishikawa and E. Ott, Chaos , 033107 (2014).[25] P. Ritchie and J. Sieber, Chaos , 093116 (2016).[26] M. I. Freidlin and A. D. Wentzell, Random Perturbations of Dynamical Systems, 2nd ed.(Springer-Verlag, New York, 1998) p. 430.[27] A. D. Ventcel’, Teor. Verojatnost. i Primenen. , 235 (1976).[28] R. Graham and T. T´el, Phys. Rev. Lett. , 9 (1984).[29] R. Graham and F. Moss, Theory of Continuous Fokker–Planck Systems , 225 (1989).[30] M. I. Dykman and M. A. Krivoglaz, Physica A , 480 (1980).[31] M. I. Dykman, M. M. Millonas, and V. N. Smelyanskiy, Phys. Lett. A , 53 (1994).
32] M. I. Dykman, I. B. Schwartz, and A. S. Landsman, Phys. Rev. Lett. , 078101 (2008).[33] H. R. Jauslin, Physica A , 179 (1987).[34] R. S. Maier and D. L. Stein, Phys. Rev. Lett. , 1783 (1993).[35] R. S. Maier and D. L. Stein, SIAM J. Appl. Math. , 752 (1997).[36] H. Touchette, Phys. Rep. , 1 (2009).[37] A. Kamenev, Field theory of non-equilibrium systems (Cambridge University Press, 2011).[38] L. S. Tsimring and A. Pikovsky, Phys. Rev. Lett. , 250602 (2001).[39] C. Masoller, Phys. Rev. Lett. , 020601 (2003).[40] B. Lindley, L. Mier-y Teran-Romero, and I. B. Schwartz, inAmerican Control Conference (ACC), 2013 (IEEE, 2013) pp. 4587–4591.[41] I. B. Schwartz, L. Billings, T. W. Carr, and M. I. Dykman, Phys. Rev. E , 012139 (2015).[42] Y. N. Kyrychko, K. B. Blyuss, and E. Sch¨oll, Eur. Phys. J. B Condens. Matter Phys. , 307(2011).[43] Y. N. Kyrychko, K. B. Blyuss, and E. Sch¨oll, Phil. Trans. R. Soc. A , 20120466 (2013).[44] Y. N. Kyrychko, K. B. Blyuss, and E. Sch¨oll, Chaos , 043117 (2014).[45] R. Sipahi, F. Atay, and S.-I. Niculescu, SIAM J. Appl.Math. , 738 (2008).[46] A. Campbell and R. Jessop, Math. Model. Nat. Phenom. , 1 (2009).[47] K. Blyuss and Y. Kyrychko, Bull. Math. Biol. , 490 (2010).[48] T. Faria and S. Trofimchuk, Nonlinearity , 2457 (2010).[49] S. Gourley and J.-H. So, Proc. R. Soc. Edinb. , 527 (2003).[50] R. P. Feynman and A. R. Hibbs, Quantum Mechanics and Path Integrals (McGraw–Hill,1965).[51] C. W. Eurich, A. Thiel, and L. Fahse, Phys. Rev. Lett. , 158104 (2005).[52] P. Kl¨oden and E. Platen, Numerical solutions of stochastic differential equations (Springer,New York, 1992).[53] S. Iacus, Simulation and Inference for Stochastic Differential Equations (Springer, New York,2008).[54] B. Abdo, E. Segev, O. Shtempluck, and E. Buks, J. Appl. Phys. , 083909 (2007)., 083909 (2007).