Abstract
In his seminal 1952 paper, Alan Turing predicted that diffusion could spontaneously drive an initially uniform solution of reacting chemicals to develop stable spatially periodic concentration patterns. It took nearly 40 years before the first two unquestionable experimental demonstrations of such reaction–diffusion patterns could be made in isothermal single phase reaction systems. The number of these examples stagnated for nearly 20 years. We recently proposed a design method that made their number increase to six in less than 3 years. In this report, we formally justify our original semiempirical method and support the approach with numerical simulations based on a simple but realistic kinetic model. To retain a number of basic properties of real spatial reactors but keep calculations to a minimal complexity, we introduce a new way to collapse the confined spatial direction of these reactors. Contrary to similar reduced descriptions, we take into account the effect of the geometric size in the confinement direction and the influence of the differences in the diffusion coefficient on exchange rates of species with their feed environment. We experimentally support the method by the observation of stationary patterns in redox reactions not based on oxihalogen chemistry. Emphasis is also brought on how one of these new systems can process different initial conditions and memorize them in the form of localized patterns of different geometries.
1. Introduction
In ‘The Chemical Basis of Morphogenesis’ Alan Turing theoretically predicts that the coupling of reaction processes with the diffusive transport of species can spontaneously lead to the development of spatially ordered patterns. In Turing's own words ‘Such a (reaction–diffusion) system, although it may originally be quite homogeneous, may later develop a pattern or structure owing to an instability of the homogeneous equilibrium, which is triggered off by random disturbances’ [1]. In this investigation, mainly focused on the onset of instabilities, he shows that these can take six essentially different forms. As he himself pointed out, the most interesting form is the one which leads to ‘stationary waves’. This investigation, proposed as a simplified explanation for pattern development in living systems, is a milestone in both theoretical biology and applied mathematics. Once developed, the concepts have received an evergrowing interest from many disciplines over the last 50 years. His unique published contribution to the field of reaction–diffusion (RD) processes is his most cited article with more than 3600 citations. Nowadays, it receives about 200 citations per year (figure 1) [2]. Interestingly, in the early years, his article attracted more attention from biologists than from applied mathematicians, the citations being essentially found in biological journals (embryology, botany and genetics) [2]. The real attraction to his work came at the beginning of the 1970s, as seen by the sudden increase of citation rate in figure 1, as his approach received theoretical thermodynamic support from Prigogine and his coworkers [3–5]. They clearly pointed out that the direct spontaneous transition from a uniform state to a stationary patterned state (the Turing bifurcation) requires that the system involves at least one positive (e.g. autocatalysis) and one negative feedback (inhibition) process, includes chemical species with appropriately different diffusion coefficients and, last but not least, operates far from thermodynamic equilibrium. The above theoretical developments came nearly simultaneously with the first published observations by Zhabotinsky et al. [6,7], and later by Winfree [8,9], of longlasting periodic wave trains in batch solutions of the now socalled BZ reaction [6,10].
Turing's approach to pattern formation was theoretically further extended in the 1970s and 1980s to pattern developments far beyond the onset of a spatial instability in an initially uniform system. These early developments were particularly made popular by the work of Hans Meinhardt [11] and James Murray [12]. Theoretical investigation of RD processes is still an advanced field in applied mathematics. Besides nonlinear chemical [13,14] and biochemical systems [11,12,15,16], this field now embraces domains as diverse as low pressure flames [17], electronhole plasmas in semiconductors [18], gas discharge devices [19], surface catalysis [20], electrochemistry [21], nonlinear optics, [22,23], defect distribution in materials irradiated with energetic particles or light [24], spreading of vegetation [25,26].
The indisputable demonstration of Turing's predictions in chemical systems was provided by the ‘Non Linear Dynamic Systems’ group in Bordeaux, at the end of 1989 (figure 2), nearly 40 years after their theoretical prediction [27]. This breakthrough was made possible by operating a newly discovered reaction—the chloriteiodidemalonicacid (CIMA) reaction [28]—in newly invented open spatial reactors of different geometries. Soon after the observation of stationary hexagonal and striped pattern arrays in the CIMA reaction [29], dramatic static and spatiotemporal patterns were observed in the ferrocyanideiodatesulphite (FIS) reaction [30], an other newly discovered bistable and oscillatory reaction [31]. In this case, the spacefilling stationary patterns resulted from interacting fronts induced by finite amplitude perturbations. These commonly form labyrinthine structures, as illustrated in figure 3. Besides labyrinths, travelling and birth–death dynamics of ‘spots’ were observed in this system [32,33]. The combined complementary observations considerably boosted theoretical investigations on RD processes [34–40]. By the 1990s, hundreds of wellcharacterized oscillatory reactions [13,14] which, in principle, have suitable kinetics for producing nontrivial spatial patterns, were known but the above two reactions remained the only examples of singlephase isothermal aqueous solution systems to produce stationary RD patterns for nearly 20 years.
Why such a standstill? We can put forward two reasons:

— the quantity of induced new experiments kept several groups busy for more than a decade. These included the studies of bifurcation sequences between different twodimensional and threedimensional pattern planforms [29,41,42], the effect of parameter ramps on pattern selection and orientation [41,43,44], (figure 4 illustrates a pattern bifurcation unfolding in a ramp), pattern modes resulting from the interaction between a spatial Turing and different uniform instabilities [45–47], the role of the chemical nature of the supporting gel [46–49], the external control of Turing pattern wavelength and geometry [50,51], the effect of moving boundaries on pattern selections [52], etc.;

— though the initial discoveries resulted from targeted research, they benefited from favourable circumstances. Full understanding of the origin of the necessary appropriate differences in the diffusion coefficients [53,54] and the actual complexity of open spatial reactors was gradually understood [55,56]. Only then could an effective design method be proposed to discover stationary patterns in different reactions [57].
In this paper, we briefly introduce the theoretical background on which our semiempirical method for the design of stationary patterns relies. The method is first exemplified by numerical simulations based on a simple but realistic model proposed to describe a class of oscillatory chemical reactions. To retain a number of important features of real open spatial reactors, but keep calculations to a minimal complexity in the numerical illustrations, we propose a collapsed dimension RD formulation that improves a similar commonly used CFUR approximation [58]. Then, we provide experimental support for the method by presenting the design of stationary patterns in the first two halogenfree redox reaction systems. In a final part, we use the peculiar properties of some localized stripe structures, developed by one of the systems, to produce geometrically controlled localized filament patterns. A number of technical details are reported in the electronic supplementary material.
2. The basis of an experimental design method to produce reaction–diffusion patterns
For the sake of mathematical tractability and minimal dynamic variable requisite, theoretical studies on the selforganizing properties of RD systems are commonly developed on a set of two coupled partial differential equations of the following form (2.1) and (2.2). 2.1and 2.2where u is an activator (∂_{u}f_{u} > 0)) and v is an inhibitor (), and are the characteristic length scales and and are the characteristic timescales of the activator and the inhibitor, respectively. The properties of the patterns and pattern formation scenarios in these systems are primarily determined by the parameters and . Importantly, it has been shown that, whatever the route to the formation of stationary patterns is (spontaneous or resulting from finite amplitude localized perturbations) in systems with steadyuniform boundary conditions, their developments are always promoted when [34,40]. These conclusions extend to systems with more than two variables but then, what has to be considered is the relative spread and timescales of the activatory process rather than that of a particular selfactivatory species. Any species interfering with this process could do the job.
A general description of an open real chemical RD system involves the local reaction terms and the linear diffusive (Fickian) transport of the chemicals across the system. The necessary matter exchanges between the system and the environment are described by the boundary conditions. The chemical kinetic term is usually a complicated function that contains dozens of reaction steps and rate constants. Even if all the information were known, which is very rarely the case, it would be quite tedious to find the appropriate patterning regions in such a large set of parameters, since these systems are usually not analytically solvable. However, if we assume that the dynamics is essentially of activatory–inhibitory type, then, according to theoretical works [34,40], we anticipate that we need only to control the main timescales, including the matter exchange process, and the length scales of the diffusive transport. The appropriate timescale of matter exchange must be in the range of that driving the reactions, since if the system is dominated by the exchange rate, the concentrations and spatial state properties do not differ from that of the environment. On the opposite limit, in case of slow exchange compared with the timescale of the activatory reaction, the system approaches thermodynamic equilibrium. In both limits, the development of temporal and/or spatial symmetry breaking patterns are excluded. Here, it is shown that we only need three carefully chosen experimental parameters, to control the relative values of the time and length scales, and to locate the suitable conditions for pattern formation.
Our design method starts with the selection of an autoactivated reaction supplemented by an appropriate antagonist chemical feedback which would lead, when operated in a spatially homogeneous continuous stirred tank reactor (CSTR), to a steadystate bistability parameter domain and to one of large amplitude oscillations. In an ideal CSTR, the dynamics of such systems can be described by: 2.3where describes the local kinetics, k_{0} is the rate constant of the matter exchange, 's and 's are the concentrations in the reactor and in the input flow, respectively.
In pertinent sections of the phasediagram, the domains of bistability and oscillations should exchange at the crossing of the stability limits of two steady states. Typically, one of the steady states has a low extent of reaction () and is referred to as the state, while the other steady state, referred to as the state, is characterized by a high extent of reaction, (the 's refer to the concentrations at thermodynamic equilibrium in a closed system). As illustrated in figure 5, such ‘crossshape phasediagrams’ [59] enable step 1, a systemic identification of parameters and , which most suitably act on the timescale of the activatory process, and on the amplitude and timescale of the negative feedback process. As a function of appropriate feed species, all known oscillatory chemical reactions can be made to exhibit such type of phasediagram.
In step 2, the spatially distributed RD dynamics of the selected reaction is explored likewise in an unstirred onesidefedreactor (OSFR). Such open spatial reactors commonly consist in a thin disc of gel in contact with one face with the steady uniform content of a CSTR. It is a RD system strongly confined in the direction orthogonal to the feed surface and extended in the two others (for details see §§1 & 6 in the electronic supplementary material). In this reactor configuration, the feed composition at the gel/CSTR interface is fixed (Dirichlet boundary condition) by the CSTR composition, while noflux boundary conditions apply at the opposite surface of the gel. The CSTR composition is described by equation (2.3) if the feedback of the gel contents on the state of the CSTR is neglected. While the concentration field within the gel is described by the following equation: 2.4
The state of the gel content is indirectly controlled by the CSTR feed composition, . If the CSTR is on the state, the gel is fed with a ‘reacted’ solution of chemicals and no further significant chemical changes are expected within the gel (, ). We refer to this state as the Tstate of the gel content. In contrast, different stationary and oscillatory states can develop in the gel if the CSTR content belongs to the state branch. The gel is thereby fed with a quasiunreacted steady composition.
Quite naturally, OSFR studies restrict to the domain of feed parameters, where this state branch is stable. All other CSTR state branches are excluded from this work. In this condition, two different steadystate profiles, not breaking the boundary symmetry, may develop orthogonally to the feed surface in the OSFR. Typically, one is characterized by a relatively low extent of reaction along the confinement direction. In this Fstate of the gel, the OSFR composition differs everywhere only little from that of the CSTR (, ). The other steadystate profile is characterized by a low extent of reaction at the contact with the CSTR content and a large extent of reaction on the opposite side. This state, which exhibits a more or less sharp switch between a low and a high extent of reaction composition, is referred to as the mixed or Mstate of the gel [56] (for illustration see §1 of the electronic supplementary material). Spatial bistability corresponds to the range of feed parameters for which the stability domain of the F and Mstates overlap. In addition to the above steady profiles, spatiotemporal oscillations or waves can be also be observed in the gel part, when in contact with this stationary ℱstate.
We anticipate that the OSFR phasediagram in the parameter subspace of the above RD systems (equation (2.4)) will qualitatively follow that of the CSTR (equation (2.3)) if the timescale of matter exchanges of the principal chemicals within the OSFR are of the order of , the residence time of the CSTR. This is used as a rule of thumb in the experimental work.
The presence of longrange inhibitory and shortrange activatory processes, owing to species with appropriate differences of diffusion coefficients, considerably enhance the capacity of RD systems to develop stationary RD patterns. If such diffusive transport differences do not preexist in the system, they have to be induced by physical or chemical additives. This step 3 is a most challenging experimental problem to solve. A direct way is to use a physically modified medium where the relative diffusive transports of species are modified. This was successfully obtained in a BZ reaction system when the homogeneous aqueous medium was replaced by a microemulsion, where nanosize aqueous droplets are dispersed in a continuum oil phase. Polar species, such as the autoactivatory bromous acid, remain trapped in the droplets and present a much smaller effective diffusion than the apolar ones, like bromine a precursor of bromide ions the main inhibitory species, which are able to diffuse across the oil phase [60]. This success was made possible by the exceptionally longlasting batch oscillatory properties of this reaction. A more widely applicable chemical method in singlephase chemical systems is to introduce a low mobility (e.g. high molecular weight) complexing agent that selectively and reversibly binds the activator. It is theoretically well established that such types of reversible binding leads to an apparent slowing down of diffusivity of the activator, the lengthening of the timescale of the activatory process [53,61] and also, but usually neglected, to a lengthening of the exchange time of the activator with the environment.
The effect of increasing amounts of complexing agent () is primarily studied in the region of and , where spatiotemporal oscillations are initially observed. This choice is guided by the theoretical proof that, by adding a low mobility complexing agent, a Turing bifurcation can be reached only in the region of phasespace, where the initial uniform steady state is an unstable node [61]. As is increased, the extent of the region with periodic oscillations decreases, more complex spatiotemporal behaviours emerge and, above a critical value, the condition for the spontaneous development of stationary patterns is reached. Stationary patterns can also be obtained outside the initially oscillatory domain, but in this case, finite amplitude local (temporary or permanent) perturbations are required to initiate the development of localized or nontrivial space filling patterns.
3. Numerical illustration of the design method
In this section, we use an extended version of a general model proposed by Rábai [62] for a family of approximatively 20 variants of pH oscillators. The parameters of the model were chosen in chemically realistic ranges of values but were not adjusted to account for specific experimental observations. The calculations should be considered as wellcontrolled experiments on an ideal RD system showing how our design method works and giving access to information that is very difficult to obtain directly in real laboratory experiments (for technical details see §4 in the electronic supplementary material). The model consists of four reaction steps (R1)–(R4). R1 R2 R3 R4
It is based on the acidcatalysed oxidation of a weak acid (HA) by an oxidant (B) in reaction (R 2), supplemented with a protonconsuming reaction (R 3) induced by component C. Component S^{−} is a macromolecule that binds the activator (H^{+}) in reaction (R 4), and has significantly lower diffusivity than all the other species. P and Q are end products. The corresponding set of differential equations can be found in the electronic supplementary material. We use experimentally manageable parameters to explore the RD system: , and are direct functions of the input CSTR feed concentrations of HA, C and S^{−}, respectively.
In step 1, a typical CSTR crossshape phasediagram is obtained in the parameter section as shown in figure 5a. At low , the system reaches a steady ℱstate branch characterized by a high pH (low acid concentration). At high , the system reaches the steady 𝒯state branch characterized by a low pH (high activator concentration). When is low, the range of stability of these two steadystate branches overlap in the bistable region. Figure 5b illustrates the pH steadystate values of these branches in the bistability region. and are the critical control parameter values for which the ℱstate and 𝒯state, respectively, loose stability. Kinetic oscillations require a negative feedback process with a large enough amplitude and suitable timescale separation. This is reached for . A typical pH versus time curve obtained in these conditions is shown in figure 5c.
To simulate step 2, we use a new simplified description of an OSFR. The confinement direction x is collapsed but, contrary to a commonly used similar approximation [58], we introduce an effective matter exchange rate with the environment (CSTR) that takes into account the mobility differences of the species and the effect of L_{x} (the thickness). Arguments for this ‘differentiatedfeedunstirredreactor’ (DFUR) approximation are discussed in §2 of the electronic supplementary material. In this heuristic approximation, the dynamics of a twodimensional (x, y) gel strip OSFR is thereby reduced to onedimensional calculations described by equation (3.1): 3.1
The confinement direction x is squeezed out, thus the RD dynamics of the gel content is calculated only along yaxis.
For comparison with the CSTR uniformstate phase diagram, we first calculate the uniform state or zerodimensional phase diagram of equation (3.1). The size of the system in the feed direction (the thickness, L_{x}) plays a crucial role in the capacity of the system to develop nontrivial patterns. For an evaluation of appropriate values, the competition between the exchange time of the activator inside the gel and the timescale of the activatory process is considered: This leads to . If we take , we obtain with our selected values of parameters (see §3 in the electronic supplementary material) that L_{x} should be of the order of 0.07 cm. This value falls in the range of gel thickness used in the experiments.
Figure 6 shows a corresponding uniform state OSFR phase diagram for L_{x} = 0.1 cm. At low χ_{1}, the content of the OSFR differs little from that of the CSTR (): the extent of reaction is low and corresponds to a high pH Fstate branch (figure 6) of the gel. At high χ_{1}, the extent of reaction is high at x = L_{x}, e.g. the proton concentration is about three orders of magnitude higher than at the contact with the CSTR. this low pH state is referred as the M state of the gel (figure 6). For , the stability domains of these two states overlap, the system is spatially bistable. Note the topological similarities between the CSTR and OSFR phasediagrams. Typical pH values of the steadystate OSFR spatial bistability branches are represented in figure 6b. For oscillations develop in the OSFR part. Figure 6c provides a typical example of pH oscillations obtained with our model approximation.
The model system is based on an activatory species that has a much higher diffusion coefficient than any other species in the system and thus is inappropriate, as is, for producing stationary RD patterns. We proceed to step 3 and explore the effect of adding a low mobility complexing agent S for the activator. Typical sets of uniformstate phase diagram sections (χ_{1}, χ_{3}), at fixed values of χ_{2}, are shown in figures 7 and 8 for the CSTR and OSFR contents, respectively. Figures 7a and 8a illustrate the effect of S^{−} in the bistability domains (). The extent of the χ_{1} bistability range initially increases with χ_{3} for both the CSTR and the OSFR. Then, the and values in the OSFR tend to become independent of χ_{3}, while the CSTR bistability limits slowly shift to higher χ_{1} values. More significant changes are observed when the complexing agent is introduced in the parameter region where oscillations are initially observed (figures 7b and 8b). The CSTR oscillations vanish above and bistability resumes. This infers that above the necessary timescale separation between the activatory and the inhibitory processes for the oscillations is no longer fulfilled [59]. A similar phenomenon is observed in the gel, but the critical is higher than . The changes are essentially owing to the shift to lower χ_{1} values of the stability limits of the Tstate and Mstate of the CSTR and OSFR, respectively.
We now test the capacity of our RD model system for producing stationary and nonstationary onedimensional patterns at a few space points in selected regions of the uniformstate phase diagram. In the spatial bistability domain () and in the absence of the complexing agent (χ = 0), appropriate finite perturbations can trigger an interface between the F and Mstates. The direction of propagation of this interface can be controlled by χ_{1}. At high χ_{1}, the Mstate invades the Fstate, we say we have a (+)front. Conversely, at low χ_{1}, the Fstate is dominant, we then say we have a (−)front. Whatever the type of front, they annihilate in headon collisions (figure 9a). No stationary pattern is found in these conditions.
If χ_{3} > 0, three different scenarios can be distinguished:

— at χ_{2} = 0 (no inhibitory process), the above triggered fronts continue to annihilate in headon collisions (figure 9b). Even if they are considerably slower, they were never found to be able to form a stable pair in these conditions. However, stable pairs of fronts and other patterns are obtained in the presence of supercritical amounts of both the inhibitor C~ χ_{2} and the complexing agent S^{−}~ χ_{3}: in the range of spatial bistability, close to the lower stability limit of the Mstate (χ_{1} = 0.05), triggered (−)fronts moving in opposite directions can form a stationary localized pulse as shown in figure 9c. The probe for stable front pairing is a major phase of this last stage of the method [63], in the absence of spontaneous pattern developments;

— at () and (), large amplitude relaxation spatiotemporal oscillations are typically observed (figure 10). These oscillations typically show brief uniform periodic switches to low pH (~4) and return back to a high pH (~7.5) plateau;

— at and , the spontaneous development of stationary patterns is observed. Figure 11a shows such a transition from the homogeneous Fstate when χ_{1} is changed from 0.12 to 0.13. First the gel content switches to the homogeneous Mstate, but this state is spatially unstable and a stationary pattern spontaneously develops after a while (figure 11a). This transition is subcritical: the back change is obtained when χ_{1} is changed from 0.11 to 0.10. The same type of stationary pattern is reached starting from the homogeneous Mstate at and decreasing χ_{1} to 0.15 (figure 11b). This transition is also subcritical: the back transition occurs when χ_{1} is increased from 0.17 to 0.18.
Figure 12 gathers the computed stationary spatial concentration profiles of the active species of the model. It shows that though species B and C play an important kinetic role in the system, their relative spatial amplitude modulations are about 60 times smaller than that of [A_{tot,G}] = [A^{−}] + [AH], the other reactant driving the activatory process (figure 12b). This points to a possibility for reducing the number of dynamic variables in the R_{1}–R_{4} model by setting B and C constant, for this set of kinetic parameters. Other important information is that the concentration of free protons H^{+} is everywhere much lower than that of the acid forms HA and HS (figure 12a). Notably, they are predominantly trapped in the very low mobility HS species only at the peak concentration of H^{+} owing to the sharp drop of [A_{tot,G}].
While the observation of stationary pattern is the target of the method, the above model also shows many other complex spatiotemporal patterns that will be discussed in a forthcoming publication.
4. Experimental illustrations
Using the abovedescribed semiempirical method, we have recently managed to experimentally observe stationary RD patterns for the first time in two halogen free redox reactions [64]. These reactions are based on the acid autocatalytic oxidation of sulphite to sulphate by hydrogen peroxide (HP) as described by the following chemical steps: R5and R6where the rate of step (R 6) has been shown to take the following form: [62]. It is a close approximation of step (R 2), where HSO would correspond to AH. In a CSTR, this HPS reaction exhibits bistability between an alkaline (pH ~8–9) unreacted state branch and an acidic (pH ~4–5) reacted state branch. In an OSFR, this system exhibits spatial bistability between an Fstate and an Mstate (figure 13a).
In a disc OSFR, under the spatial bistability conditions, an interface between the low and high acid steady states can be initiated by a local perturbation or by the particular conditions found at the rim of the disc reactor (see §5 in the electronic supplementary material). Figure 13b shows a time–space plot of the collision between two (−)fronts initiated at the rim of the disc OSFR. As in theoretical calculations (figure 9a), no stationary or spatiotemporal symmetry breaking pattern of any kind was observed in the absence of a low mobility complexing agent.
In accordance with the above general scheme, this autocatalytic reaction must be supplemented by a negative feedback reaction to produce stationary patterns. In a first example, a proton sink is provided by the oxidation reaction of ferrocyanide to ferricyanide by HP: R7
This step is a variant of (R 3), where the oxidant (B in R 2) is involved in the negative feedback process. The input feed concentration of ferrocyanide plays the role of parameter χ_{2} in the model. Sections of the CSTR and OSFR phase diagrams of this HPSF reaction are shown in figure 13a. As required, the CSTR and OSFR bistability domains shrink with increasing . Above a critical level, the CSTR contents exhibits large pH oscillations, similar to those reported by Rábai [62]. Above a slightly higher critical value, the OSFR exhibits spatiotemporal oscillations (figure 13c). The diagrams are in qualitative agreement with those obtained with model (figure 6).
At this stage, sodium polyacrylate (NaPAA) with molecular weight of 15 000 D is added to the system (). This weak acid with pK_{a} values typically ranging between 6 and 5 will reversibly capture protons in the reacted state: R8
This low mobility weak acid was shown to be effective for producing stationary RD patterns in other acid autoactivated reactions [54,57]. Figure 14 illustrates the expected development of stationary patterns in the presence of a supercritical [NaPAA]_{0} feed for the same value of other feed parameters as in figure 13c.
Hydrogen carbonate provides a nonredox alternative reaction to remove hydrogen ions through the equilibrium reactions (R 9) and (R 10): R9and R10
Note that in this HPSC reaction, the sum of steps (R 9) and (R 10) resumes to (R 3), here [HCO_{3}^{−}]. In this reaction, the timescale of the negative feedback is independent from the concentration of HP [65]. The remarkable difference from the previous case is that on gradually introducing hydrogen carbonate, the onset of an oscillatory dynamics is immediately accompanied by complex spatiotemporal patterns.
Starting from the uniform Fstate and increasing at fixed low values of induces the development of complex spatiotemporal patterns mixing travelling domains and quasistationary filamentous structures of Mstate. The latter are permanently erased by the passing of a travelling domain but reconstruct at a different position behind. The general aspects of this spatiotemporal pattern are illustrated in figure 15a. On increasing , the size of travelling domains decreases while the density of the filamentous structure increases (figure 15c). The space–time plots taken along a vertical median in the previous figures show that these spatiotemporal dynamics are quite ‘chaotic’. The addition of mM can lead to expected stationary patterns [64]. This is less than half the amount needed to stabilize stationary patterns in the HPSF system.
According to the general theoretical analysis, the observation of symmetry breaking spatiotemporal patterns in the HPSC system in the absence of NaPAA suggest that the space scale of the activatory process is already significantly shorter than the antagonist processes. In this system, the effective production of protons during one kinetic cycle is lower than in any other documented pH patterning reaction. Thus, the contribution of other lowmobility protonbinding species such as the colour indicator could have an important effect on the HPSC spatiotemporal dynamics (see electronic supplementary material in Szalai et al. [64]).
It was also found that the dynamics of this system is very sensitive on the thickness of the gel, L_{x}. A 10 per cent gradual difference between two parts of the disc can be dramatic. We make use of this sensitivity to unfold in space a sequence of pattern bifurcations in this system, as in figure 4. Using this fact, figure 16 reveals that doubling the BTB indicator concentration can indeed lead to stationary cellular patterns, in the absence of NaPAA. The thickness of the disc is still close to 0.5 mm everywhere, but the top part is very slightly thicker (~10%) than the bottom part. From the bottom, there is first a transition from an Fstate to complex spatiotemporal patterns mixing waves and filamentous structures as in figure 15. As we go up, the density of filaments grows and a transition to a stationary cellular pattern is reached. The latter spread over a narrow band in the middle of disc. Then, separated by a very narrow region of spatiotemporal patterns reminiscent of a TuringHopf mixed mode [37], a domain of travelling lowamplitude phasewaves is observed before reaching the stable uniform Mstate at the top of the disc. Such spatial unfolding of bifurcation should help to understand the as yet theoretically undocumented spatiotemporal patterns presented in figure 15.
5. Pattern development beyond the Turing instability
In recent years, RD studies have developed beyond the understanding of pattern morphogenesis. They are proposed to solve logic [66] and pathfinding problems [67,68] or ways to process and store information [69]. In this last context, localized RD structures are playing a privileged role. Contrarily to standard Turing patterns that are usually spacefilling patterns with a relatively welldefined wavelength, localized patterns are isolated structures of one state embedded in a ‘sea’ of another state. Such structures are usually found in bistable systems. In the following, we report on the development of isolated low pH stripe or filament structures in a uniform high pH background.
In the HPSF reaction, at not too high values of but high value of [NaPAA], there is no oscillation and no spontaneous development of pattern. In these conditions, starting from the uniform Mstate a decrease of [H_{2}SO_{4}]_{0} ~ χ_{1} destabilizes this state below a critical value. Interestingly, this destabilization of the Mstate always occurs first at the rim of the disc and a growing annular domain of Fstate slowly invades the disc and eventually fills it. No localized Mstate pattern (e.g. isolated spot) is produced. But, if in these feed conditions, two disconnected domains of Fstate grow, their merging can be impeded. The (−)fronts can be made to repel each other sufficiently to form stable filaments. A switch from the M to the Fstate can be induced by a light pulse at wavelengths corresponding to the Fe(CN) absorption band [70]. We use this property to produce and study the stability of filament patterns induced by localized perturbations of different geometries (for technical details, see the electronic supplementary material).
All the following experiments start with similar initial conditions: the whole disc is set uniformly in the Mstate and the acid feed is gradually decreased until a spontaneous switch to the Fstate starts at the rim of the disc (see §5 in the electronic supplementary material). Then, different geometric localized light perturbations are made, ‘the initial event’. This initial event is ‘processed’ and ‘memorized’ in a stable pattern. In the simplest case, the light perturbation is made at the centre of the shrinking Mstate domain, as shown in figure 17a. Subsequently, a domain of Fstate grows at the centre and comes to a close encounter with the external Fstate annulus. The two (−)fronts stop when they are at a distance of about 2 mm and form a relatively stable circular Mstate filament. The space–time plot in figure 17b shows this sudden stop of the colliding (−)fronts which initially propagated at 0.04 mm min^{−1}. If the acid feed is further decreased by 2 per cent, the Mstate pattern collapses quasiuniformly and we are left with a uniform Fstate. Conversely, if the acid feed is increased by 2 per cent, the ring pattern swells, a pair of (+)fronts develop and the Mstate reinvests the whole disc, the memory is reinitialized. A new initial perturbation can be tested after decreasing back the [H_{2}SO_{4}]_{0} value.
Figure 18 shows the situation when four approximatively equally spaced point perturbations forming the vertices of a square are made. The system first produces a fourleaf clover pattern of Mstate filament with a ‘crossing’ of nearly 90° at the centre. This central quadrangular crossing is unstable and, on a longer timescale, it splits in two ternary connections that gradually separate and start to equalize. The angles between the filaments then slowly approach 120°. Figure 18b,c show the space–time plots taken along a vertical and horizontal line approximatively through the centre of the disc.
In figure 19a, the initial geometric perturbation consists of two horizontal lines and two spots near the ends of the central Mstate band. The filaments with free tips formed at the top and bottom parts shrink by retraction of the tips (retraction rate mm min^{−1}). While the central long filament is stabilized by the two end loops. The time–space plot taken horizontally just above the final central filament shows (figure 19b) that the filament junctions with the loops move apart at a rate of 0.015 mm min^{−1}. The vertical diameter of the end loops shrink at mm min^{−1}, i.e. the loops slowly and uniformly shrink.
The behaviour of these large cellular structures bear strong similarities to surface tensiondriven structures such as in soap foams. However, the very longterm stability of these patterns could not be accurately studied owing to a slow hydrolysis of the ferrocyanide stock solution that induces a very slow increase of the pH in the CSTR. We may wonder if analogue types of information processing and memory mechanism could be at work in the brain of living systems.
Note that the present filament structures show no transversal nor elongation instabilities—i.e. do not break into rows of spots, develop no meanders, neither side nor tip budding phenomena. These particular properties, which have not been encountered in other solutionphase chemical patterning systems as yet, can also be used for solving an ‘entrance–exit’ problem in a maze. A simple such example can be found in §6 of the electronic supplementary material.
6. Concluding remarks
Even if our semiempirical design method is successful [64,71], a more fundamental theoretical analysis is required to better set its limits. In this perspective, a more detailed numerical exploration of the patterning capacity of the pH oscillator model (R1)–(R4) will be pursued. The use of our DFUR reduced dimensional description of real OSFRs, which better preserves the actual physics of these open spatial reactors than the similar commonly used CFUR approximation [58], will reliably help these investigations. It should make possible to account for a broader range of actual experimental observations and effectively guide the exploration of the patterning capacities of still more reactions. Why is it interesting to multiply the number of experimental examples? Contrary to formal theoretical and numerical studies where all the parameters of a model are expandable, so that all the possible patterning capacities of a model can be explored, real chemical reactions have physical constraints. Only a small part of a potential phase diagram can be reached by one reaction. Different reactions allow access to different combinations of parameter values. New patterns and new routes to patterns can be discovered, especially since real experimental systems are multitime and multispace scale systems. RD systems with large numbers of variables are still seldom explored in theoretical studies. The observations of different routes to stationary patterns and other symmetry breaking spatiotemporal patterns in a variety of inorganic reactions should also pave the way for the observation of such patterns in biologically meaningful reactions and thus further test Turing's original idea on biological morphogenesis.
Acknowledgments
We acknowledge supports from the CNRS and the Hungarian Research Fund OTKA (77986, 100891). D.C. and N.T. acknowledge supports from the ESF research networking programme ‘FUNCDYN’. J.H. is supported by the KTIAOTKAEU 7FP ‘Mobility’ programme (80244) and N.T. is supported by the Hungarian Development Bank. We also thank P. Borckmans at the ULB, Brussels, for fruitful discussions.
Footnotes
One contribution of 13 to a Theme Issue ‘Computability and the Turing centenary’.
 Received January 16, 2012.
 Accepted March 6, 2012.
 This journal is © 2012 The Royal Society