Cite as: Phy. Rev. E (2025); https://doi.org/10.1103/v1h2-42w1
PDF
version
Dynamical systems with time delays have garnered significant attention across various scientific and technological domains due to their intricate behaviours and inherent complexity . Time delays, which capture the non-instantaneous nature of interactions, have been integrated into mathematical models in diverse fields, including biology , chemistry , optics , economics , and ecology . They are especially relevant in contexts such as enzyme reactions , cell signalling , and coupled oscillatory reactions, where temporal delays critically shape system dynamics. These systems, which closely resemble real-world scenarios, have become pivotal in both theoretical investigations and practical applications, particularly in generating chaotic signals and advancing the development of chaotic systems. Time delays profoundly influence system dynamics, driving phenomena such as chaos in one-dimensional systems , amplitude death , and oscillation death . Moreover, they can give rise to complex behaviours, including chaotic synchronization and chimera states in coupled systems, as well as spatio-temporal pattern formation and the triggering of extreme events . These diverse and intricate phenomena underscore the critical role of time delays in shaping the rich dynamics of interacting systems.
In chemical kinetics, certain reactions may require a finite time interval before influencing subsequent reactions. When this delay becomes significant, it must be incorporated into the mathematical model of chemical reactions using delay differential equations. Such time delays, often manifested as delay feedback, are prevalent in chemical and biological systems and can impact both reaction and diffusion processes within the reaction-diffusion model, which have been extensively studied, especially in the context of pattern formation, target pattern formation , spontaneous motion , spatiotemporal instabilities , breathing oscillations , and oscillatory behaviours. In the study of chemical models, the Brusselator is a well-known theoretical framework for oscillatory chemical reactions that has been studied under the influence of diffusion and time-delayed feedback, with studies including Hopf bifurcation analysis , pattern formation , feedback control , and the emergence of extreme events .
In addition, extreme events have garnered significant attention, leading to extensive exploration across various dynamical systems and providing crucial insights into the mechanisms driving these phenomena . Hence, examining extreme events within the context of time-delayed systems is essential to understanding their occurrence and potential impacts in natural settings. The study of extreme events has been observed in various time delay systems, including a neuronal model coupled with two time delays , and has been shown to be significantly affected by increasing delay feedback strength, which can increase the occurrence of events in an optical system , while it reduces and eventually eliminates the occurrence of extreme events in a forced oscillatory system , highlighting the relevance of delay effects in their dynamics.
Motivated by these observations, this paper focuses on extreme events in the Brusselator model, specifically investigating their occurrence under the influence of time delay feedback in the absence of diffusion, and further extending the study to the presence of diffusion to account for spatial interactions . In addition, we identified transient chaos phenomena in the proposed system. To validate these numerical findings, we implemented an experimental circuit as an alternative approach to explore and confirm the practical applicability of the theoretical dynamical model. The circuit dynamics were analysed in the laboratory and compared against numerical simulations of a normalized model. This normalized model was developed through appropriate substitutions of circuit variables and parameter transformations .
This paper is structured as follows: In Sec. [Model Description], we introduce the mathematical model and perform a linear stability analysis. Sec. [Numerical Results], focuses on the numerical investigation of the system, demonstrating the occurrence of extreme events through time series analysis, probability distribution functions, and return period plots. The characteristics of transient chaos are examined in Sec. [Transient chaos], and the numerically obtained results are validated by experimental results conducted using an analog electronic circuit in Sec. 4. The role of diffusion in spatial interactions is discussed in Sec.[Influence of Diffusion]. Finally, we summarise our findings in Sec. [conclusion].
An autocatalytic chemical reaction is considered within the two-dimensional theoretical Brusselator model, incorporating self-feedback of the delayed autocatalyst concentration . The Brusselator reaction’s schematic diagram , shown in Fig. 1. We have highlighted the time-delayed autocatalyst self-feedback loop in red. This model corresponds to the reactive scheme (Fig. 1) which, after suitable rescaling, results in dimensionless equations. The theoretical Brusselator model governing equations with self-feedback delay is expressed as follows, \[\begin{aligned} \nonumber \dot{x} &=& a - b x - x + x^2 y + \epsilon x_{\tau}\\ \dot{y} &=& b x - x^2 y , \label{equation_1}\end{aligned}\] where \(x\) and \(y\) are the state variables of the system, representing the dimensionless concentrations of autocatalysts, while \(a\) and \(b\) are constant system parameters corresponding to the dimensionless concentrations of reactants.
The delayed variable \(x_\tau\) denotes the delayed concentration of the autocatalyst \(x(t - \tau)\), with \(\tau\) as the time delay and \(\epsilon\) as the delay feedback strength in the reaction.
The equilibrium point of the system described by Eq. ([equation_1]) is determined by constant concentrations of autocatalysts over time, obtained by setting the time derivatives of \(x\) and \(y\) equal to zero and at steady state, where \(x(t - \tau) = x(t)\). As a result, the system possesses unique equilibrium points given by \((x^*, y^*)~=~(\frac{a}{1-\epsilon}, ~\frac{b(1-\epsilon)}{a})\). The stability of a stable equilibrium point (\(x^*, y^*\)) is typically defined as a condition in which all nearby trajectories asymptotically converge to this point \(t \to \infty\), and hence it is analysed by introducing a small perturbation around the equilibrium point (\(x^*, y^*\)) with a time-dependent function.
The stability analysis is based on the characteristic equation, with eigenvalues determined from \(|J_0~+~e^{-\lambda\tau}~J_\tau~-~\lambda\text{I}|~=~0\), where \(J_0\) represents the Jacobian with respect to the equilibrium point (\(x^*, y^*\)), and \(J_\tau\) is the Jacobian with respect to delay (\(x_{\tau}\)). The characteristic equation derived for the equilibrium point can be expressed as: \[\lambda^2 + \lambda \left( \frac{a^2}{(1 - \epsilon)^2} - b + 1 - \epsilon e^{-\lambda \tau} \right) + \frac{a^2}{(1 - \epsilon)^2} \left( 1 - \epsilon e^{-\lambda \tau} \right) = 0.\] This characteristic equation, expressed in transcendental form, was analysed using the Newton-Raphson method to numerically determine the minimum of twenty roots (eigenvalues). If all eigenvalues have negative real parts, the equilibrium point is \(\it{stable}\); however, if any eigenvalue has a positive real part, the equilibrium is \(\it{unstable}\).
The local stability of the system as a function of delay time (\(\tau\)) is depicted in Fig. 2(a) by the maximum real part of the eigenvalue (\(\Re(\lambda_{max})\)), showing that the system transitions from an unstable equilibrium to a stable one through a Hopf bifurcation point with stable equilibrium points depicted in orange, unstable points in blue, and the Hopf bifurcation point is highlighted in red, indicating the critical transition between the stability regimes. Similarly, in Fig. 2(b), the maximum real part of the eigenvalue (\(\Re(\lambda_{max})\)) is depicted as a function of the feedback strength \(\epsilon\), showing that the system remains primarily in a stable regime but exhibits intermittent unstable points for \(\epsilon \in [0.4975,0.5839]\), with transitions into and out of unstable regimes through Hopf bifurcation points.
The stability plot for the two-parameter (\(\tau,\epsilon\)) plane, corresponding to the eigenvalues \(\lambda\), is depicted in Fig. 3. The system parameters are fixed at \(a=0.21\) and \(b=0.686\), while delay time (\(\tau\)) is varied from 0 to 0.2, and the delay feedback strength (\(\epsilon\)) is adjusted from 0 to 1. Hence, the system exhibits a confined unstable regime at intermediate delay feedback strength, with the Hopf bifurcation curve serving as the boundary between stable and unstable regions. Different colour zones in Fig. 3 indicate various stability regimes, such as the stable regime shown in orange, the unstable regime indicated by blue, and the Hopf bifurcation curve highlighted in red.
The investigation begins with a bifurcation analysis, examining the dynamic response of the system ([equation_1]) as the time delay is varied within the interval \(\tau \in [0, 0.2]\). For the numerical analysis, we employ the fifth-order Runge-Kutta method through the method of steps in the Julia programming language with a fixed step size of 0.01, using constant system parameters \(a = 0.21\), \(b = 0.5\), \(\epsilon = 0.5\) and time delay (\(\tau\)) is varied to study the system dynamics. The initial condition is set to, \((x_{0},~y_{0}) = (0.01, ~0.02)\) and the initial history function is defined to be the same as the initial condition for \(\forall t \in [-\tau, 0]\).
The one-parameter bifurcation diagram is presented in Fig. 4, showing the system response for the autocatalyst concentration (\(x_{max}\)) as a function of the time-delay parameter (\(\tau\)). Notably, significant amplification in the system amplitude is observed at \(\tau = 0.085\), indicating the emergence of extreme events. This abrupt transition denotes the onset of rare, large-amplitude excursions within the system. These extreme events persist across a range of delays until \(\tau = 0.115\), beyond which, the system transitions out of this extreme event regime to bounded chaos.
To quantify extreme events deviations from nominal chaos, several statistical approaches are utilised, including the (\(90^{th}-99^{th}\)) percentile of the probability distribution of events or using a threshold derived from event magnitude, such as the peak-over threshold method . Here, we calculated the threshold height as \(H = \langle x_{\text{peak}} \rangle \pm m~\sigma_{x_{\text{peak}}}\), where \(\langle x_{\text{peak}} \rangle\) represents the mean of the peak state variable, \(\sigma_{x_{\text{peak}}}\) is the standard deviation of \(x_{\text{peak}}\), in this context, the term peak refers to the local maxima and minima of the state variable, and \(m\) is a multiplicative integer specific to the system, determining how far an event deviates from the mean value, we set \(m = \pm 8\) to qualify events as extreme, as depicted in Fig. 4, where these events are highlighted in green. Extreme events are further distinguished by the red in the bifurcation diagram, emphasising their significance in the system’s dynamics.
To explore the occurrence of extreme events, the numerically obtained time series of the state variable \(x(t)\) and its corresponding phase portrait are shown in Fig. [time_series]. These illustrate both non-extreme and extreme events scenarios derived from the one-parameter bifurcation analysis as a function of time delay, with all other system parameters kept consistent with the previously utilised values. For a time delay \(\tau=0.123\), the system exhibits the bounded chaos within the system’s state space, which is depicted in Fig. [time_series]a(i), while the corresponding time series is presented in Fig. [time_series]a(ii) which exhibits bounded chaos and lacks the ability to surpass the threshold line(\(H\)) for both the local maxima and minima of the state variable, depicted in the horizontal green and red lines, respectively. When the time delay parameter \(\tau\) is varied to \(0.1\), as illustrated in Fig. [time_series]b(i), the system exhibits sudden, large excursions in the state space. The trajectory deviates from the bounded chaotic regime, undergoing a significant excursion before eventually returning to the system’s bounded chaotic attractor. This rare and occasional large deviation is classified as an extreme events (EEs). The corresponding time series, shown in Fig. [time_series]b(ii), reflects the deviation as large amplitude excursions that exceed the threshold line in both the local maxima and the minima of the state variable, characterising them as extreme events.
In addition, a projection of the 3D phase plot is depicted in Fig. [Statistical characterisation] to illustrate how the bounded chaotic attractor in Fig. [Statistical characterisation]a(i) expands into a larger exclusion zone as extreme events in Fig. [Statistical characterisation]b(i), by including the time delay of the state variable \(x_{\tau}\) as the third axis.
The probability distribution function (PDF) of the state variable \(x\) was analysed to confirm the infrequent occurrence of extreme events (EEs). Figure [Statistical characterisation](ii) illustrates the PDF calculated over an extended time span of \(3 \times 10^8\) time units, ensuring the system has evolved beyond transient states. For non-extreme events, the PDF remains confined within a bounded region below the threshold line, as shown in Fig. [Statistical characterisation]a(ii). Conversely, for extreme events, the distributions exhibit extended tails with values exceeding the threshold (\(H\)), marked by a dashed vertical line in Fig. [Statistical characterisation]b(ii). This highlights the presence of extreme events and demonstrates their low probability of occurrence beyond the defined threshold.
In Fig. [Statistical characterisation](iii), depicts the return period (RP) for both non-extreme and extreme events, with the threshold \(H\) indicated as a dashed horizontal line. The return period quantifies the average interval between successive occurrences of an event and is calculated using \(RP = \frac{1}{1 - p(x_{\text{peak}})}\), where \(p(x_{\text{peak}})\) represents the empirical cumulative distribution function, which defines the probability of a value being less than or equal to \(x_{\text{peak}}\). The figure highlights that non-extreme events do not exceed the threshold (\(H\)) in Fig. [Statistical characterisation]a(iii), while extreme events surpassing \(H\) exhibit notably high return periods, ranging from \(2075\) to \(109979\) in Fig. [Statistical characterisation]b(iii). These values correspond to exceptionally low probabilities of occurrence (\(1 / RP\)) between \(0.0048\) and \(0.9 \times 10^{-5}\). This analysis emphasises the rarity and exceptional nature of extreme events, confirming their rare occurrence within the system.
The distinct relationship between the occurrence of extreme events and the delay feedback strength is revealed in Fig. 5 which shows the number of extreme events, defined as peaks surpassing the threshold \(H\), plotted as a function of the feedback strength \(\epsilon\). Specifically, for both weak and strong delay feedback strengths, the system does not exhibit any extreme events. However, as the feedback strength reaches an intermediate range, the number of extreme events sharply increases, indicating that extreme events are most likely to occur within a specific range \(\in [0.39, 0.58]\), of moderate feedback strengths. This underscores the critical role that feedback strength plays in modulating the emergence of extreme events, highlighting that extreme behaviours are confined to a particular delay feedback regime.
Further, the system’s overall dynamics across the parameter space (\(\tau,\epsilon\)) are illustrated in Fig. 7(a), with \(\tau\) ranging from 0 to 0.5 and \(\epsilon\) spanning from 0 to 1. Our objective is to identify and differentiate regions where the extreme events (EEs) occur from those with non-extreme events (NEEs). To distinguish the EEs regions as a function of the parameters \(\tau\) and \(\epsilon\), we determined a threshold height (\(H\)). During a long simulation of \(10^8\) time units, if the maximum values of the \(x\)-variable (\(x_{\max}\)) exceed the threshold height for \(m = 8\), the region is categorized as experiencing EEs (marked by red points in Fig. 7(a)). However, if the system remains below \(H\), those parameter values are classified as part of the NEE region (depicted by grey points in Fig. 7(a)). Through this analysis, we confirmed the presence of EEs within the specified ranges of \(\tau\) and \(\epsilon\) used time series data. In Fig. 7(b), where the data are interpolated for clarity, it is evident that for this range \(\epsilon \in [0.35, 0.7]\), there is a significant probability (red-orange colour region) of extreme events occurring in this range. From Fig. 5, we conclude that EEs occur only within a particular range of moderate feedback strengths.
In this section, the phenomenon of transient chaos, characterized by the appearance of chaos with a finite transient lifetime , as observed in the delay system governed by Eq. [equation_1] retained the same set of parameter values while systematically varying the time delay precisely to \(\tau = 0.38855\). The system undergoes an extended transient state characterized by chaotic behaviour, persisting for approximately \(t \approx 1.9 \times 10^4\) iterations. Beyond this state, the system transitions into an asymptotic state marked by stable 4-T periodic cycles.
Fig. [tc]a(i) illustrates the phase portrait of the system during transient chaos, where the chaotic dynamics are depicted in blue, and the eventual periodic state is represented in red. The temporal evolution of the state variable \(x\) is shown in Fig. [tc]a(ii), capturing the chaotic behaviour during the transient state. The subplot emphasises the 4-T periodic state observed after the transition. The geometric characteristics of the chaotic transient phase are further explored on the Poincaré surface of section in the (\(x\), \(y\)) plane, as presented in Fig. [tc]a(iii). During the transient behaviour, the chaotic dynamics are represented by a collection of randomly distributed points in blue, while the periodic (4-T) state is indicated by four distinct points in red.
Additionally, the temporal Poincaré plot for the state variable \(x\) is depicted in Fig. [tc]a(iii) and a(iv), providing further insights into the system’s behaviour. The chaotic transient dynamics dominate the plot initially, transitioning into a clearly discernible 4-T periodic state, which is highlighted in the subplot. This result underscores the intricate dynamics of transient chaos and the system’s eventual stabilisation into a periodic state while varying time delay.
This section describes the design of a simple analog circuit implemented in the laboratory to validate the numerical results of the two-dimensional non-autonomous time-delayed Brusselator model outlined in Eq. [equation_1]. The schematic of the experimental Brusselator circuit with self-feedback delay is presented in Fig. [circuit diagram]. For a comprehensive discussion of the implementation and a detailed description of the Brusselator oscillator circuit, readers are referred to Ref. . We have fixed the values of other circuit elements as follows: The capacitance values of the capacitors are fixed at \(C_1\), \(C_2 =\)2.2 nF. The resistances of resistors are set as follows \(R_1\), \(R_3\), \(R_6\), \(R_7\) and \(R_8=\) 10 k\(\Omega\), \(R_2=\)6.8 k\(\Omega\), \(R_5=\)3.3 k\(\Omega\), \(R_4 =\) 100 k\(\Omega\). The circuit is powered by a dual-voltage power supply of \(\pm\)12 V.
In Fig. [circuit diagram], the time-delay unit, highlighted within the red dashed box, is implemented using a Sallen-Key low-pass active filter configured as a low-pass Bessel filter to achieve the desired delay . The time-delay block adopts an op-amp-based approach, following the methodology described by Buscarino \(et~ al\). , employing a cascade of Bessel filters. This design is particularly effective for generating time delays in the millisecond range, as required for this study.
The time-delay unit comprises a cascade of two second-order low-pass Bessel filters, each constructed using the Sallen-Key topology. The discrete components utilised in the circuit are resistors \(R_9\), \(R_{10}\), \(R_{11}\), \(R_{12} = 10\) k\(\Omega\) and capacitors \(C_3\), \(C_4\), \(C_5\) \(C_6 = 10\) nF. Each filter introduces a measured time delay of 1.3 ms, achieved by appropriately tuning the resistor (\(R_9\) to \(R_{12}\)) and capacitor (\(C_3\) to \(C_6\)) values. This corresponds to a dimensionless delay (\(\tau\)), calculated as: \[\begin{aligned} \tau = \frac{T_{delay}}{RC}=1.3. \nonumber\end{aligned}\] The described circuit was constructed using readily available discrete components. The time-delay unit was designed with two cascaded filters, enabling the effective time delay to be adjusted. This configuration allowed for a systematic analysis of the circuit’s behaviour across varying delay values.
With the nominal system parameters of the Sallen-Key blocks (specifically employing two Sallen-Key stages), the circuit exhibited the anticipated bounded chaotic behaviour and extreme events. Experimental results are presented in Fig. [Experimental observations of EEs], which shows a digital oscilloscope trace of the normal chaotic attractor (Fig. [Experimental observations of EEs](a)) and extreme events (Fig. [Experimental observations of EEs](b)) in the phase plane formed by the voltages across capacitors \(C_1\) and \(C_2\) (i.e., the \(v_{c_1}\)–\(v_{c_2}\) plane), where \(v_{c_1}\) and \(v_{c_2}\) correspond to the variables \(x\) and \(y\), respectively, in the numerical simulations. The corresponding time-delay waveform in the (\(v_{C_2}(t-\tau\))) plane is shown in Fig. [Experimental observations of EEs](c). The experimental observations align closely with the numerical simulations, validating the model and circuit design.
In this section, we extend our analysis by incorporating diffusion into the system, thereby introducing spatial interactions that significantly influence and modify the system’s dynamical behaviours. In particular, we examine how diffusion impact the emergence of extreme events and compare these results with the delay-induced extreme events observed in the non-diffusive case.
With the addition of diffusion, the Brusselator model extends to a general reaction-diffusion form , expressed as \[\begin{aligned} \nonumber \dot{x} &=& D_X \nabla^2 x + a - b x - x + x^2 y + \epsilon x_{\tau}\\ \dot{y} &=& D_Y \nabla^2 y + b x - x^2 y , \label{equation_diffusion}\end{aligned}\] where \(D_X\) and \(D_Y\) are the diffusion coefficients for the \(x\) and \(y\) respectively, while \(\nabla^2 x\) and \(\nabla^2 y\) represent their spatial diffusion.
The initial condition of system ([equation_diffusion]) is set by introducing small random perturbations around the homogeneous steady state as \[x_{(X, Y, 0)} = x^* + \kappa \eta_{(X,Y)}, ~y_{(X, Y, 0)} = y^* + \kappa \eta_{(X,Y)},\] where \(x^* = a/(1-\epsilon)\) and \(y^* = b(1-\epsilon)/a\) are the homogeneous steady-state solutions, \(\kappa\) is a small perturbation parameter (\(\kappa = 0.01\)), and \(\eta_{(X, Y)}\) is a spatially varying random function sampled from a uniform distribution over \([0,1]\).
The system parameters are taken to be the same as in Sec. [Numerical Results], with diffusion coefficients set as \(D_X=1\) and \(D_Y=1\). The spatial domain is discretized using a uniform grid of \(N_X \times N_Y = 20 \times 20\) points, with spatial step sizes \(\Delta X = \Delta Y = 1\). The system ([equation_diffusion]) is subject to Neumann boundary conditions, \[\frac{\partial X}{\partial n} = 0 ~,~ \frac{\partial Y}{\partial n} = 0 ~\text{on}~ \partial \Omega,\] ensuring zero flux across the boundaries of the spatial domain.
For the diffusion case as well, prominent high-amplitude extreme events are observed in the time series of \(x_{(X_{10}, Y_{10}, t)}\), which exceed the predefined threshold (\(H_{max}~\And~H_{mini}\)) for \(m = 4\) and are therefore classified as extreme events, as shown in Fig. [Diffusion](a). The spatiotemporal evolution of \(x_{(X, Y, t)}\) along the fixed spatial slice \(Y = 10\) during the extreme event interval highlighted in gray in Fig. [Diffusion](a)—is depicted in Fig. [Diffusion](b).
These results suggest that while time delay alone is sufficient to induce extreme events in the non-diffusive system, the interplay between delay and diffusion does not significantly alter the occurrence of extreme events.
In this study, we explored the dynamics of the theoretical Brusselator chemical model with the incorporation of self-feedback of autocatalyst (\(x\)) through a time delay (\(x(t - \tau)\)), uncovering several intriguing phenomena, including extreme events and transient chaos. Both time delay (\(\tau\)) and delay feedback strength (\(\epsilon\)) were shown to play a critical role in the emergence of extreme events (EEs), as confirmed through statistical measures. Our analysis revealed that the EEs predominantly occur within specific regions of the system parameter space. By systematically varying the time delay, we observed transitions in the system dynamics from bounded chaotic behaviour to the emergence of EEs and transient chaos. The steady state of the system was found to be highly dependent on the interplay between the system parameters (\(a \And b\)) and the delay feedback strength. Notably, varying the feedback strength altered the attractors in phase space and introduced instability in intermediate feedback regions, where the number of EEs was observed to be noticeably high. These findings were corroborated using a two-parameter phase diagram, highlighting the prominent regions associated with the occurrence of EEs. Parallely, the real-time experimental circuit is constructed for the self-feedback delayed system and obtains similar results from numerical observations. In the experiments, we have created the variable delay and obtained the chaos and extreme event oscillations. By incorporating diffusion into the model, we analyzed the spatial evolution of extreme events and assessed the role of spatial interactions, and we carried out several numerical and statistical measures to confirm the obtained phenomena.
The proposed model, especially with delay, provides a simple yet powerful platform to explore the emergence of complex behaviours such as chaos and multistability. It bridges fields like physics, biology, and engineering, demonstrating universal principles of nonlinear and delayed systems. Feedback delay models are instrumental in designing and understanding control systems, particularly in contexts where delays are unavoidable (e.g., communication networks, robotics). The modified Brusselator with delay could contribute to the development of digital twins for dynamic systems with time-lagged interactions. It is a simple, analytically tractable model that allows researchers to test new mathematical methods for analysing nonlinear and delayed systems.
The work of S.V.M was catalyzed and financially supported by the Tamil Nadu State Council for Science and Technology, Department of Higher Education, Government of Tamil Nadu, India. K.T. acknowledges the Chennai Institute of Technology for the experimental and computational assistance. SS acknowledges the Basic Research Program of the National Research University, Higher School of Economics, Moscow.
Data available on request from the authors.