Stationary localized structures and the effect of the delayed feedback in the Brusselator model

1. Introduction

The spontaneous emergence of dissipative structures in chemical systems arises from a principle of self-organization that can be either in space and/or in time [1–3]. Prigogine demonstrated that chemical oscillations are perfectly compatible with the second principle of thermodynamics extended to far from equilibrium systems (see overviews on this issue [4,5]). Spatial dissipative structures generally evolve on macroscopic scales and can only be maintained by continuous application of a non-equilibrium constraint. A classic example of spatial self-organization has been provided by Turing [6]. This mathematical prediction has been supported by Prigogine's explanation of the physical mechanisms underlying the formation of dissipative structures. The Turing–Prigogine dissipative structures are characterized by an intrinsic wavelength. The formation of spatial structures is attributed to the balance between a nonlinear process originating from a chemical reaction that tends to amplify spatial fluctuations and a diffusion process that tends, on the contrary, to restore uniformity. When a system operates far from equilibrium, dissipation plays an important role in the stabilization of self-organized structures. Dissipative structures in chemical systems have been experimentally observed thanks to the development of chemical open reactors with a Chlorite–Iodide–Malonic-Acid reaction [7,8]. The above mechanism has been widely used to explain the formation of stationary periodic dissipative structures not only in chemistry but in various biological and physical systems [9–20].

Besides the spatially periodic structures, the same mechanism predicts the possible existence of aperiodic and localized structures. The latter consist of isolated or randomly distributed spots surrounded by regions of the uniform state [21–25]. They correspond to the confinement of energy, chemical concentration, light or biomass density in one or more spatial dimensions (see recent overviews on this issue [15,16,26,27]). A vast amount of theoretical work has been carried out on localized structures in chemical reaction–diffusion systems [28–38], while there is still a lack of experimental evidence of stable localized spots. This is because isolated localized structures (LSs) often exhibit a self-replication phenomenon that affects the circular shape of the LS [39,40] and leads to the formation of an extended patterned state.

In this paper, we consider the paradigmatic Brusselator reaction–diffusion model. We report on stationary two-dimensional solutions of the model consisting of localized spots. We provide a bifurcation analysis of the localized solutions, which arise in a subcritical pitchfork bifurcation and undergo a self-replication instability leading to bound states of localized solutions. We use two-dimensional continuation techniques provided by the matlab package pde2path [41] to construct the bifurcation diagram associated with two-dimensional localized structures. In the last part of the manuscript, we investigate the effect of delayed feedback on the stability of localized structures. We show that both periodic and localized structures exhibit a spontaneous motion in an arbitrary direction, whereas in the absence of time-delayed feedback, a single spatial spot and periodic patterns remain stationary. The bifurcation diagram associated with a two-dimensional single localized spot is constructed. We characterize the motion by estimating the speed of moving dissipative structures and by analysing the threshold associated with the onset of motion as a function of the feedback parameters.

The paper is organized as follows: after the introduction, we analyse the bifurcation structure of localized spots in §2. In §3, we implement delayed feedback control. We analyse the impact of delayed feedback on the linear regime. In §4, we characterize the transition from a stationary to a moving localized spot. In §5, we conclude.

2. Localized structures in the Brusselator model

The Brusselator model is a classical reaction–diffusion system, proposed by Prigogine & Lefever [2]. The evolution equations of the Brusselator model read

Here, X = X(x, y) and Y = Y (x, y) denote the concentrations of the interacting chemical species, whereas A and B are externally controlled concentrations which play the role of control parameters. The parameter D is the ratio between the diffusion coefficients D_y and D_x of the Y and X concentrations, respectively.

The homogeneous steady state of (2.1) is (X_s, Y_s) = (A, B/A). It can undergo a Turing–Prigogine bifurcation leading to the formation of a stationary oscillation in space or an Andronov–Hopf instability leading to time-periodic chemical oscillations. The threshold of the Turing–Prigogine bifurcation is given by , D⁻¹ < 1 whereas the Andronov–Hopf instability occurs if B > 1 + A². In what follows, we focus on the case of stationary spatially localized patterns, so that we choose B < 1 + A² so that the system operates far from the Andronov–Hopf bifurcation for the stability of (X_s, Y_s).

The formation of spatially periodic structures such as hexagons or stripes and transitions between them have been intensively investigated for the model (2.1) [42]. Besides periodic patterns, the same mechanism predicts the possible existence of aperiodic, localized structures. They consist of isolated or randomly distributed localized spots embedded in the uniform state [34,37,38]. Under uniform chemical constraints, non-equilibrium reaction–diffusion systems which exhibit a subcritical Turing–Prigogine instability display stable localized spot patterns in the pinning parameter range where there is coexistence between two different states: a uniform branch of stationary state solutions and a branch of spatially periodic solutions. In this hysteresis loop, a single circular spot consisting of a pinning front between a stable uniform solution and, e.g. a hexagonal branch of solutions, is stable [37]. Indeed, the Brusselator model equation (2.1) supports stable two-dimensional localized solutions with periodic boundaries in both x and y directions as shown in figure 1 (2). One can observe these structures in direct numerical simulations of equation (2.1); however, this tool seems unfit for a systematic analysis of the bifurcational structure. We therefore deploy numerical continuation techniques provided by the matlab package pde2path [41]. Using this package, we are able to track stable and unstable solutions as well as folds and other bifurcations. Providing this analysis in two spatial dimensions is a novelty in theoretical biology; however, it has been successfully applied in other fields of research, such as modelling of thin films [43,44] and nonlinear optics [45].

• Download figure
• Open in new tab
• Download PowerPoint

Starting from a stable localized solution obtained by direct numerical simulations with periodic boundaries, we explore the parameter space using the parameter B as our primary continuation parameter. The L₁-norm of the activator field X, i.e.

with being the mean value of the field, is used as a measure to plot bifurcation diagrams. Since the structures we are interested in possess a radial symmetry, we perform our calculations on a quarter of the real domain with Neumann boundaries for the sake of numerical performance and then expand the plotted solutions shown in figures 1 and 2. For our calculations, we use a mesh adaption method both to refine the mesh where it is required and to coarse the mesh if possible.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

The result of this bifurcation analysis is depicted in figure 3. The branch of localized solutions bifurcates from the branch of quadratic periodic solutions in a supercritical pitchfork bifurcation, very close to the critical Turing-point. Recent work on this topic [46] suggests that this is a finite size effect and that, in an infinite system, both branches would start at the critical point. At first, an unstable low amplitude solution starts to arise (cf. figure 1, first row) that is, however, still unstable. Along the branch indicated by the blue dotted line in figure 3 the structure grows in amplitude and finally reaches stability in a fold at B ≈ 0.4. This is the main solution branch we are interested in, where a stable localized solution (cf. figure 1, second row) exists approximately in the range of 0.4≤B≤0.59 (solid blue line). This regime of existence is in good agreement with results from direct numerical simulations. At B ≈ 0.59, the localized solution loses its stability in a subcritical pitchfork bifurcation. Further along the now again unstable branch (blue dotted line), the structure develops a dip at the centre and transforms into a ring-like structure (cf. figure 1, third row). It then runs through further bifurcations and becomes even more unstable, so we stop to trace the branch at B ≈ 0.71. Also, at the afore-mentioned bifurcation point, another unstable solution branches off (green dotted line). On this branch, the radial symmetry of the solution is broken. The solution shows a self-replicated structure of two bound localized structures in the first component X of our system. In the faster diffusing second component Y , the solution looks more like an elliptic localized solution (cf. figure 1, fourth row). This solution branch declines down to the periodic solution, suggesting that it connects with a periodic branch; however, up to now we are not able to determine the type of bifurcation at the end of this branch.

Aside from the single localized solution, we were able to track another branch of localized solutions, also branching off from the periodic orbit in a similar way as the single solution. The bifurcation diagram is depicted in figure 4. In this case, two equal localized spots form a bound solution with a distance of half the domain size between them. The solution starts again as an unstable low-amplitude solution (cf. figure 2, first row) and then reaches stability (solid blue line) through a fold and a nearby bifurcation, leading to the described bound state (cf. figure 2, second row). At this bifurcation, another solution branches off (dotted green line). Plotting this solution shows that on this branch one of the two spots is decreasing (figure 2, fourth row). Therefore, the discrete translation symmetry of half the domain length is broken for solutions on this emerging branch. Similar to the single spot solution, the bound solution loses its stability in a second bifurcation and leads to a ring-like solution (figure 2, third row) which becomes more unstable through further bifurcations (not shown). At the second bifurcation point, another solution, where each spot is replicated, branches off (dotted green line). The solution is depicted in the fifth row of figure 2 and refers to the breaking of the radial symmetry. Although we were also able to find this stable localized solution of two spots in direct numerical simulations, one should note that this solution depends on the domain size. Contrary to the single localized solution, the exact shape changes with varying domain size since, due to the periodic boundaries, every solution possesses two neighbours that interact with each other.

All described solution branches are depicted in figure 5. One could suspect that the decreasing spot solution (figure 2, fourth row) branching off near the fold of the bounded localized solution connects to the single spot solution. But this is not the case since the decreasing spot does not vanish completely. Instead, it remains a bound solution of a smaller and a larger spot. One observes that the corresponding solution branch runs similarly through further bifurcations as the single localized spot solution does. However, in contrast to the single localized solution, it is unstable in the whole parameter space.

• Download figure
• Open in new tab
• Download PowerPoint

This section shows that, using numerical continuation techniques in two spatial dimensions, we are able to analyse the bifurcation structure of localized solutions. Most important for further investigations of localized solutions in the Brusselator model are the two regimes where stable localized solutions exist.

3. Branches of stationary solutions for the Brusselator modelwith delayed feedback

In another line of research, it has been shown that time-delayed feedback can strongly impact the space–time dynamics of both periodic and localized states in various systems out of equilibrium described by the Swift–Hohenberg equation with time-delayed feedback via the so-called Pyragas control [47–53]. In chemical reactions without diffusion, delay cannot be ignored in industrial applications where recycling of unreacted reagents reduces the cost of the reaction. The output stream of a continuously stirred tank reactor is sent through a separation process [54,55]. In reaction with diffusion and with global delayed feedback, a photoemission electron microscope was used to continuously image lateral distributions of adsorbed species, and feedback is introduced by making the instantaneous dosing rate of the catalytic carbon monoxide oxidation dependent on real-time properties of the imaged concentration patterns [56].

The effect of spatial inhomogeneities for the delayed Swift–Hohenberg equation has been investigated both analytically and numerically in [57] and a transition from oscillating to depinning solutions has been characterized. In this case, the time-delayed feedback acts as a driving force [57]. The effect of noise responsible for the formation of dissipative structures in the delayed Swift–Hohenberg equation has been investigated recently in depth by Kuske et al. [58]. Furthermore, it was shown that delayed feedback can induce motion and breathing localized structures in chemical reaction–diffusion systems [52,59]. In advanced photonic devices, the effect of the phase on the self-mobility of dissipative localized structures has been theoretically investigated in [60–62]. In driven Kerr cavities described by the Lugiato–Lefever equation, time-delayed feedback induces a drift of localized structures, and the route to spatiotemporal chaos has been discussed in [63–66].

Recently, time-delayed feedback control has attracted a lot of interest in various fields of nonlinear science such as nonlinear optics, fibre optics, biology, ecology, fluid mechanics, granular matter, plant ecology (see recent overview [67]), and the excellent book by Erneux [54]. Numerical simulations of the Brusselator model with delayed feedback using Pyragas control have provided evidence of moving periodic structures in the form of stripes and hexagons [68] or superlattices [69]. Dissipative structures and coherence resonance in the stochastic Swift–Hohenberg equation with Pyragas control have also been investigated [70].

We incorporate in the Brusselator model (2.1) a time-delayed feedback term. The resulting evolution equations of the Brusselator model with time-delayed feedback read [37]

The delayed feedback parameters are the strengths of the feedbacks (η_x, η_y) and the delay time τ. The approximation we use to model the delay term is valid in the limit where the strength of the delay is small. In this limit, the feedback concentrations are small enough for their evolution to be described by a single delay term.

The homogeneous steady states with delayed feedback obey

The homogeneous steady states may be bistable. In what follows we focus our analysis on the monostable regime. The linear stability analysis of the homogeneous steady states with respect to perturbations of the form leads to the transcendental characteristic equation First, let us examine the case of a Turing–Prigogine instability where the eigenvalue σ vanishes for a finite wavenumber q = q_c. Above the threshold associated with this instability, there exists a finite band of Fourier modes q − ² < q² < q + ² with where F = D(B + 1 − η_x − 2X_sY_s) + X²_s − η_y. The Fourier modes are linearly unstable and trigger the spontaneous evolution of the state variable towards a stationary, spatially periodic distribution of the chemical concentration, which occupies the whole space available in the chemical open reactor. The critical wavenumber at the onset of the Turing–Prigogine instability reads and the threshold is a solution of the following equation The marginal curves associated with the Turing–Prigogine instability are shown in figure 6 for different values of the feedback strength, η_x = − η_y (the instability region is on the right-hand side of the curves). For negative values of η_x, the delay destabilizes the system and enables the Turing–Prigogine instability to occur for B < B_TP = (1 + A/D)², where B_TP is the threshold in the absence of delayed feedback. The threshold is thus shifted towards smaller B. For positive values of η_x, the system restabilizes itself and only a finite domain of B provides unstable Fourier modes. Increasing the delay strength further leads to the disappearance of the unstable band.

• Download figure
• Open in new tab
• Download PowerPoint

In figure 7a, we present a one-dimensional bifurcation diagram for max(X) as a function of the feedback strength η for A = 0.6, B = 0.8, D = 150 and feedback parameters η_x = η_y = η and τ = 10. Here, the green line represents the branch of localized solutions with a spatial profile as shown in figure 7b. Numerical continuation of this branch, as a function of the feedback strength, has been carried out by the Newton method that allows finding both stable (represented by solid line) and unstable (represented by dashed line) states. The branch of localized solutions is limited by an Andronov–Hopf bifurcation (diamond symbol) and saddle-node (square symbol) bifurcations. The localized solution emerges from the homogeneous steady state represented by the blue line (solid and dashed for the stable and unstable parts, respectively), which is delimited by a subcritical Turing–Prigogine bifurcation (star symbol). Characteristic of this system is that it is highly multistable: not only a region of simultaneous stability of the localized and the homogeneous solution exists in figure 7a, but there also exist stable branches of spatially periodic solutions with different wavelengths. Examples of two branches are shown in figure 7a represented by the red and the magenta lines. They are obtained by a direct numerical integration of equations (3.1). The corresponding spatial profiles of X are shown in figure 7c. Note that there exist several branches of moving LSs consisting of more peaks that are not shown in figure 7.

• Download figure
• Open in new tab
• Download PowerPoint

4. Transition from stationary to moving dissipative localized structureswith delayed feedback

The delayed feedback may also induce a drift bifurcation of the localized structure in the Brusselator model, similar to that reported in the case of Swift–Hohenberg [48], Lugiato-Lefever [66] equations and the semiconductor laser model [49,60–62]. Such a feedback-induced drift of a localized structure is illustrated in the one-dimensional case in figure 8a for the case of A = 0.6, B = 4.5 and D = 150 and feedback with a fixed strength η_x =  − η_y = 0.05 and delay time τ = 200. We also illustrate the two-dimensional case. In figure 9, we present an example of drifting localized structure for Brusselator parameters of A = 0.6, B = 0.45 and D = 150 and feedback with a fixed strength η_x =  − η_y =  − 0.17 and delay time τ = 20. We consider an initial condition which consists of stationary localized structures. During time evolution, while the structure remains quite symmetric and stationary (figure 9a), it develops an asymmetric shape and starts moving along an arbitrary direction since the Brusselator model is isotropic in space as shown in the snapshots of figure 9b,c. This drift then goes on indefinitely without change in speed or direction as the LS goes through the boundary condition. Similarly, figure 10 presents an example of drifting honeycomb pattern for Brusselator parameters of A = 4.5, B = 30 and D = 8 and feedback with a fixed strength η_x =  − η_y =  − 0.1 while keeping the delay time the same as in the previous case. All numerical simulations are conducted for periodic boundary conditions. To initialize the simulation, the provided initial condition is also used as the feedback for the first τ units of time.

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

• Download figure
• Open in new tab
• Download PowerPoint

The localized structure is launched at t = 0 and starts moving after a certain timespan, which decreases as the feedback strength is increased (see also [49]). Furthermore, the localized structure moves faster as η is increased. Our bifurcation analysis demonstrates that a localized state is a particular case since it coexists with stable periodic patterns with different wavelengths. Figure 8b illustrates the case of delayed induced drift for a periodic pattern. Note, however, that asymmetric LS solutions in one dimension [30,34] and in two dimensions [71] have been reported even in the absence of the delayed feedback control.

In order to calculate the threshold of the delay-induced drift-bifurcation, we proceed in a similar way as in [48,60,66]. Slightly above this threshold, we consider a localized structure moving uniformly with a constant small velocity v = |v| and expand the model variables X and Y in power series of v: X = X₀(ξ) + v[X₁(ξ) + vX₂(ξ) + v²X₃(ξ) + · s] and the same expression for Y . Here, X₀(ξ), Y₀(ξ) is the stationary localized structure profile, ξ = r − v et, r = (x, y), and e is the unit vector in the direction of motion of the localized structure. Substituting this expansion into equation (3.1) and collecting the first-order terms in the small parameter v, we obtain

with w_x = e · ∇X₀ and w_y = e · ∇Y₀. The linear operator L is given by By applying the solvability condition to the right-hand side of equation (4.1), we get the drift instability threshold with S_x = 〈ψ^†_x, ∇X₀〉 and S_y = 〈ψ^†_y, ∇Y₀〉. Here, the eigenfunction ψ^† = (ψ^†_x, ψ^†_y)^T is the solution of the homogeneous adjoint problem L^†ψ^† = 0 and the scalar product 〈 · ,  · 〉 is defined as . For feedback strengths of η_x = η_y = η, we recover from (4.2) the threshold condition obtained earlier for the drift instability of cavity solitons in the Swift–Hohenberg equation with delayed feedback, ητ =  − 1 [48].

One can also calculate the drift threshold for the case of −η_x = η_y = η, where we obtain

with χ = 〈ψ^†_y, ∇Y₀〉/〈ψ^†_x, ∇X₀〉. In the one-dimensional case, for the parameters A = 0.6, B = 0.8 and D = 150, with a delay time of τ = 20, this expression gives an estimation for the feedback strength of η_analytical = 0.049 = η_numerical, which corresponds perfectly to the numerically estimated threshold. In the two-dimensional case, for the same parameters as figure 9, we reach the same conclusion with η _analytical = 0.050.

One can also calculate the velocity of the structure close to the bifurcation point. The derivation is very similar to the one performed [52] for the case of Pyragas control [47], so it will not be performed in full detail here. We choose the Ansatz

where q is a vector containing both components X and Y , q₀ is the stationary solution, r = (x, y)^T is a vector containing the spatial coordinates and R(t) is the position of the localized solution. Inserting this Ansatz into equation (3.1) and projecting on the adjoint eigenfunction ψ^† = (ψ^†_x, ψ^†_y)^T yields: with Expanding R(t − τ) leads to the evolution equation for the velocity:

showing the normal form of a supercritical pitchfork bifurcation with the bifurcation at ητ = − 1. In accordance with the results mentioned above, at ητ = − 1 equation (4.7) changes from the trivial stable solution v = 0 to a non-trivial solution with

i.e. the solution starts to drift. Note that this approximation is only valid close to the bifurcation point since any shape deformations due to additional stable modes were neglected in the Ansatz. This pitchfork bifurcation is confirmed by numerical simulation (figure 11).

• Download figure
• Open in new tab
• Download PowerPoint

5. Conclusion

We have investigated the formation of dissipative localized structures in the paradigmatic Brusselator model in a two-dimensional setting. Using path-continuation techniques, we have characterized the transition from a single spot towards periodic patterns through a self-replication phenomenon. The branch of localized solutions bifurcates from the quadratic periodic solutions in a supercritical pitchfork bifurcation, very close to the critical Turing–Prigogine point. The localized solution loses its stability in a subcritical pitchfork bifurcation and develops a dip at the centre, transforming it into a ring-like structure. At the pitchfork bifurcation point, another unstable solution, with broken radial symmetry, branches off. The solution shows a self-replicated structure of two bound localized structures. We have also found another branch of two localized spots with a distance of half the domain size between them. Similar to the single spot solution, this bound solution loses its stability in second bifurcations, leading to a ring-like solution on one branch and to replicated spots on the other branch.

In addition, we have investigated the influence of the delay feedback control on the dynamics of localized and periodic structures in the Brusselator model. In particular, we have shown that delay feedback can induce the motion of one- and two-dimensional localized spots in an arbitrary direction since the system is isotropic in the (x, y)-direction. When increasing further the strength of the feedback and keeping the delay time constant, the steady state chemical concentration exhibits a complex dynamical behaviour induced by the delayed feedback. In particular, it would be interesting to characterize the spatiotemporal chaos and the formation of rogue waves in two-dimensional settings. We shall come back to these questions in a separate publication.

Data accessibility

This article has no additional data.

Competing interests

We declare we have no competing interests.

Funding

F.T. received funds from the German Scholarship Foundation and the Center for Nonlinear Science Münster. M.T. received support from the Fonds National de la Recherche Scientifique (Belgium). K.P. acknowledges the Methusalem Foundation.

Footnotes

One contribution of 12 to a theme issue ‘Dissipative structures in matter out of equilibrium: from chemistry, photonics and biology (part 2)’.