Designing Stationary Reaction–Diffusion Patterns in pH Self-Activated

Article Cite This: Acc. Chem. Res. 2018, 51, 3183−3190

pubs.acs.org/accounts

Designing Stationary Reaction−Diffusion Patterns in pH SelfActivated Systems Judit Horvat́ h,*,†,§ Istvań Szalai,† and Patrick De Kepper‡ †

Institute of Chemistry, Eötvös Loránd University, P.O. Box 32, H-1518 Budapest 112, Hungary Centre de Recherche Paul Pascal, CNRS, University of Bordeaux, 115, avenue Schweitzer, F-33600 Pessac, France

‡

Acc. Chem. Res. 2018.51:3183-3190. Downloaded from pubs.acs.org by UNIV OF GOTHENBURG on 12/19/18. For personal use only.

S Supporting Information *

CONSPECTUS: Since Alan Turing’s 1952 pioneering work, reaction−diffusion (RD) processes are regarded as prototype mechanisms for pattern formation in living systems. Though suspected in many aspects of morphogenetic development, pure RD patterns have not yet been demonstrated in living organisms. The first observations of an autonomous development of stationary chemical patterns were made in the early 1990s. In this Account, we discuss the recent developments for producing stationary pH RD patterns in open spatial reactors. The theoretical analysis of the early experiments anticipated the possibility of finding Turing patterns in a wide range of oscillatory reactions if one could control the kinetic and diffusional rate of some key species. However, no experimentally effective method to produce stationary Turing patterns was attained before 2009, and the number of systems stagnated at two until then. The two precursor reaction systems benefited from unplanned favorable chemical properties of the RD media. Theoretical studies point out that appropriate diffusion rate differences are necessary to produce stationary patterns since a competition between an effective short distance self-activation and a long distance inhibitory process is required. This differential diffusion would naturally lead to differential exchange rates between the RD system and its feed environment, an aspect somewhat overlooked in theoretical and in primal experimental approaches. Our pattern design method takes this aspect into account. A slower diffusion of a self-activated species (here, protons), produced in the RD part of the spatial reactor, generates the accumulation of this species compared to the other species. This accumulation has to be at least partly compensated by an independent scavenging reaction. The above requirement naturally brought us to focus on two-substrate pH oscillatory reactions. Stationary RD patterns are now well documented in six pH driven reaction systems. Furthermore, the coupling with a pH dependent metal ion complexing agent led to stationary patterns in calcium ion concentration. Our effective semiempirical design method does not require a detailed knowledge of the reaction kinetics; thus it is applicable to a broad spectrum of reactions and even to synthetic biological systems. It is based on simple dynamic arguments and on general topological characteristics of a nonequilibrium phase diagram. We first illustrate our method with numerical simulations, based on a realistic but idealized general model of the two-substrate pH-oscillator reaction family, and provide a refined view of the topology of the resulting phase diagrams. Then, we exemplify its effectiveness by observations made in distinct pH self-activated systems. Analogies and differences between experiments and the model calculations are pointed out. Besides standard hexagonal arrays of spots and parallel stripes, hitherto undocumented dynamic phenomena, such as randomly blinking areas and complex dynamic and stationary filamentous structures, were observed. The main challenge, to find lowmobility complexing agents that would selectively and reversibly bind a species controlling the self-activatory kinetic path of the reaction, was readily overcome in multiple ways by anions of weak acids: not only by polymeric substances but in some cases by a pH color indicator or even smaller molecules, depending on their proton binding affinity.

■

INTRODUCTION

independent from the size of the reaction medium if the system is large enough. Under uniform environment, theoretical analysis shows that spontaneous symmetry-breaking patterns, whether they are stationary or nonstationary, require appropriate differences between the diffusion coefficients (Di) of some

The formation of spatially organized structures from an initially featureless state is a central problem of developmental biology and, accordingly, of the origin of life.1,2 Chemical reactions and diffusive transport of species can spontaneously result in the formation of concentration patterns as was predicted mathematically by Alan Turing.3 He pointed out that such patterns would have a “chemical wavelength” that is © 2018 American Chemical Society

Received: August 31, 2018 Published: November 9, 2018 3183

DOI: 10.1021/acs.accounts.8b00441 Acc. Chem. Res. 2018, 51, 3183−3190

Article

Accounts of Chemical Research species.4 In a two variable activator−inhibitor system, stationary reaction−diffusion (RD) patterns form if Dactivator < Dinhibitor. This is the standard way to obtain stationary Turing patterns emerging from a uniform steady state. The first observation of stationary symmetry breaking RD patterns in a chemical system was made with the CIMA reaction.5 The observed pattern planforms of high and low triiodide concentration corresponded to those predicted by Turing. Soon after these observations, stationary patterns were found in the FIS reaction, a pH driven oscillatory reaction.6−8 Though sustained stationary patterns were also discovered in more complex physicochemical media,9 the observation of such stationary patterns in simple solution chemistry remained limited to the CIMA and FIS reactions. These discoveries were the fruit of targeted research programs but benefited from peculiar properties of the systems. In the CIMA reaction, the physicochemical origin of these properties were identified by Lengyel and Epstein.10 They recognized that the appropriate differences in Dactivator and Dinhibitor in this system were mediated by the formation of a reversible, kinetically inert complex between triiodide (a species linked to the main activator, the iodide ion) and the slowly diffusing starch (initially introduced as a color indicator to view possible triiodide patterns). On the basis of this understanding, Lengyel and Epstein proposed a design method to produce stationary patterns in different chemical systems, but their original approach did not lead to the expected breakthrough. All reactions that are capable of temporal oscillations are theoretically ready to produce Turing patterns. Why not in practice? Hundreds of oscillating chemical reactions were already known. However, the large majority of the oscillatory reactions, like the FIS reaction, are operational only under continuous feeds of reactants. Very few have a sufficiently complex mechanism, like the CIMA reaction, that allows them to exhibit a finite number of oscillations under batch conditions. In an effective design method, which intends to work for a broad range of chemical reactions, the feed processes must be suitably taken into account. Nowadays, the most popular reactors used in this field are the open one-side-fed reactors (OSFRs), sketched in Figure 1 (see

mainly in the gel where the extent of reaction increases along the normal to the feed surface. The size of the gel along this direction (direction of feed), hereafter called the thickness Lx, determines the time-scale of diffusive matter exchange in an OSFR. The basic phenomenon, which appears when operating an OSFR with an autoactivated reaction, is spatial bistability: the coexistence, for the same boundary conditions, of two different spatial states in the gel, not breaking the symmetry of the boundary feed.12 This is the starting point of our survey.13 Then we show how to find spatiotemporal oscillations and then stationary RD patterns.

■

THE PRINCIPLE OF THE DESIGN METHOD To introduce the steps and to have a better understanding of the resulting phase diagrams, with a glimpse of their possible complexity, we start with the analysis of a numerical example using a realistic skeleton reaction mechanism. The details of the applied simulation methods are discussed in the Supporting Information. Step 1: Operate with Reactions That Show Bistability and Oscillations in a CSTR and Where Activatory and Inhibitory Processes Can Be Independently Tuned

We used the general model, proposed by Rábai, representing the chemistry of two-substrate pH-activated oscillatory reactions:14 H+ + A− ↔ HA

(R1)

HA + B → H+ + P

(R2)

H+ + C → Q

(R3)

H+ + S− ↔ HS

(R4)

Reaction R1 is the protonation equilibrium of a reductant A−, and reaction R2 is the oxidation of the protonated form HA of the reductant by an oxidant B. Self-activation of proton production occurs by their combination. Reaction R3 is a proton consuming reaction induced by a second substrate C. A detailed description of the chemistry and a list of the possible components was recently accounted by Epstein and coworkers.15 The proton binding agent S−, which plays a major role in the pattern design method,13 is incorporated by equilibrium R4. From a kinetic point of view, S− acts as a buffer, but due to its assumed low mobility (in the simulations we set DHS/DHA = 0.01), it also slows down the apparent diffusivity of protons, the self-activated species. Mass action kinetics is used for reactions R1, R3, and R4, while the rate law of reaction R2 reads v2 = (k2[H+] + k2′)[HA][B]. The CSTR phase diagram of reactions R1−R3 shows the typical cross-shaped topology of activatory−inhibitory systems,16 with a relatively small domain of bistability between the stationary states (F and T), which decreases and leaves the place to a growing domain of sustained oscillations as [C]0 is increased (Figure 2A).

Figure 1. Sketch of an open one-side-fed reactor. Adapted with permission from ref 11. Copyright 2015 AIP Publishing.

Step 2: Find Conditions Where Spatial Bistability and Spatiotemporal Oscillations Develop in an OSFR

details in Supporting Information) made of a thin strip or disc of relatively inert hydrogel (polyacrylamide, agarose, poly(vinyl alcohol)), which has one of its faces in contact with the chemical contents of a continuous-flow stirred tank reactor (CSTR). The CSTR serves as a reservoir of fresh reactants. This configuration allows maintaining the contents of the spatial part at controlled distance from thermodynamic equilibrium. Simultaneously, the friction in the pores inhibits convective fluid motions in the hydrogel, and molecular diffusion dominates the transport of soluble chemical species. The reaction can be made to evolve

To obtain RD phenomena in the supposed “gel” part, we keep the CSTR on the low extent of reaction stationary state, referred to as the F (flow) state. Under this condition, the RD medium is essentially fed by fresh concentrations of reactants, and two spatial states, with different concentration distributions along the feed direction, can develop. In the F-state of the gel, the extent of reaction is overall low and the composition is similar to that of the CSTR. The M (mixed) state of the gel is characterized by the outbreak of the self-activatory process, which here 3184

DOI: 10.1021/acs.accounts.8b00441 Acc. Chem. Res. 2018, 51, 3183−3190

Article

Accounts of Chemical Research

Step 3: Introduce a Low-Mobility Agent That Reversibly Binds a Species on Which the Positive-Feedback Kinetic Process Selectively Relies and Thus Slows down the Effective Diffusion of This Species

To describe the dynamics in the RD medium along a plane (y,z) perpendicular to the direction of the feed (x), we apply a simplified 2D description of an OSFR.17 Thus, the parameter values are shifted compared to the previous simulations made along the feed direction (x). Upon increase of the concentration of S− (stot, Figure 4A), the domain of oscillations reduces and vanishes, and a domain of stationary patterns develops. Note in the vicinity of the cross-point that the stability range of stationary patterns overlaps with those of the oscillations and of the 2D uniform states (F and M). These multistabilities lead to a complex dependence of pattern aspects on the history of parameter changes (Figure 4B). Starting from the uniform Fstate, as h0 is increased, first patterns with separated M-state stripes form. Further increase of h0 results in the formation of sloppy hexagonal arrays of F-state spots immersed into an Mstate background. A backward step with h0 leads to an interconnected M-state net. Patterns with separated stripes or with a net structure can develop at the same set of parameters, depending on initial conditions. Bistability also appears between the uniform M-state and the F-state spot-pattern. As the stability limit of the M-state is crossed by decreasing h0, a less interconnected net-like pattern appears. These sequences demonstrate that the actual aspect of the patterns in this range of parameters can strongly depend on the history of the system. Above a certain level of stot (e.g., 0.015) (Figure 4A), in agreement with the theoretical considerations, the stability limits of the uniform and the patterned states become independent from the value of stot.18 Analogously, the shape and wavelength of the patterns do not change significantly with stot within the stability domain. Due to the finite size and the noflux boundary conditions, a discrete spectrum of wavelengths with λ = Ly/n (where Ly is the size of the system and n is an integer number) is observed in the simulations. At parameters corresponding to Figure 4B (Ly = 12), the wavelength varies between 6 (n = 2) and 1.71 (n = 7). Simulations made on a larger domain (Ly = 26) resulted in a similar set of wavelengths, from 5.2 (n = 5) to 1.73 (n = 15), in agreement with the intrinsic character of the wavelength. Generally, λ decreases with h0. To scan the pattern forming capacity of pH self-activated systems, an appropriate set of parameters is h0 and c0 in the presence of a supercritical value of stot. A complex planform topology is obtained (Figure 4C), with overlaps between the patterned state and one or both of the uniform states or the oscillatory state: Increasing c0, as we move the system away from the cross-point, the aspect of the patterns changes from less ordered stripes or network to more ordered hexagonal arrays of spots or more parallel arrays of stripes (Figure 4D). At the same time, the wavelength of the patterns decreases. Noteworthy, at low c0, by scanning the states of the system with h0, only spatial bistability between the uniform F and M states can be observed even in the range where the patterned state can be also stable (domain d in Figure 4C). In this case, the branch of patterned state is isolated and can be reached only by starting from the monostable domain of patterns and by gradual changes of parameter values. At the value of other parameters in Figure 4C, c0 = 0.002 is the lowest stability limit for pattern development. When the domain of spatial bistability (F/M) vanishes at sufficiently high c0, a monostable domain of stationary patterns unfolds. However, further increase of c0 restores oscillations. But

Figure 2. Nonequilibrium CSTR (A) and OSFR (B) phase diagrams of the reactions R1−R3 calculated along the plane of the nondimensional input feed concentrations of H+ and C (h0,c0). Blue, F (flow) state; red, T (thermodynamic) state; yellow, oscillatory state; orange, M (mixed) state. Bicolor hatched areas are domains of bistability as indicated in the framed labels. Simulations were made along the feed direction (x) with a 1D model11 of an OSFR.

appears as a sharp even H+ concentration jump parallel to the feed boundary. These two states coexist over a finite range of parameters; this is spatial bistability. In the domain of spatial bistability, an appropriate local perturbation can trigger a propagating front between the two steady states. Depending on the feed concentration, this front can travel at the expense of one or the other state. Concomitantly, the expansion of the high H+ concentration zone in the thickness of the gel is likewise constraint sensitive (Figure 3a).

Figure 3. Simulated traveling interfaces between the two stable stationary states (F and M) (a) and time−space plot of spatiotemporal oscillations in the reaction−diffusion system under OSFR conditions (b). The grayscale corresponds to the logarithm of the nondimensional concentration of H+ (h). Adapted with permission from ref 11. Copyright 2015 AIP Publishing.

The dynamics of the RD system in the OSFR changes above a critical value of c0 (here 0.25), and spatiotemporal oscillations spontaneously develop in the RD part (Figure 3b). The corresponding nonequilibrium phase diagram (Figure 2B) shows a cross-shaped topology: the F/M bistability domain exchanges with a domain of sustained oscillations, at the crossing of their stability limits. 3185

DOI: 10.1021/acs.accounts.8b00441 Acc. Chem. Res. 2018, 51, 3183−3190

Article

Accounts of Chemical Research

Figure 4. Simulations with the nondimensional Rábai model made along the (y,z) plane with a simplified 2D description of an OSFR: (A) simulated phase diagram at c0 = 0.05; (B) states observed at stot= 0.03 at four different values of h0 (0.0130, 0.0160, 0.0185, 0.0200); (C) simulated phase diagram at stot = 0.04; (D) patterns observed at c0 = 0.20 at two different values of h0 (0.0190, 0.0210). Blue, F (flow) state; orange, M (mixed) state; yellow, spatiotemporal oscillations; purple, stationary RD patterns. Hatched areas are domains of multistability as indicated in the framed labels.

Step 2A: Finding Spatial Bistability

regardless, the pattern formation scenario starts to show the typical signs of a standard Turing bifurcation: the regular hexagonal array of spots and parallel stripes expected in a 2D system (Figure 4D).19

We start with the pH autoactivated subsystem, oxidant + sulfite, without substrate C or S−. The oxidant/sulfite ratio was set above the stoichiometric ratio to ensure the maximum production of protons through path R2 and to leave enough oxidant for C if path R3 is a redox sequence. Spatial bistability is expected at a gel thickness Lx where the diffusion (exchange) time τD is on the order of the induction time τchem of the autoactivated proton production, so that Lx = (τchemD)0.5 where D is the typical value for the diffusion coefficient of the feed species other than H+. Lx was usually 0.75 mm (Table S1) with τD ≈ 470−560 s. τCSTR was fixed, while τchem was tuned by the feed concentration of a strong acid, [H2SO4]0. The two spatial states (F and M) were monitored by the color change of an appropriate pH color indicator (Ind) fed in the solution. Strictly, this color indicator is an intervention and infiltrates an additional protonation equilibrium: Ind− + H+ ↔ HInd. Ideally, the concentration of Ind and HInd should be low, exhibit separated absorption bands (color filters in Table S1 are only transparent for HInd), and pKind should lie between the two lowest pH-values of states F and M. To get enough contrast with our cameras and filters, the indicator concentrations were not considerably lower than the proton concentrations ((1−3) × 10−4 M), a concentration comparable to [H2SO4]0 (see Table S1). The M state comprises a pH-switch somewhere in the depth of the gel (Figure 3a). Its position can be close to the back of the OSFR, so the optical path length, where HInd is the main absorbing species, can be much shorter than Lx (see Figure 3a). In reactions with iodate as oxidant, the intermediate species triiodide (iodine + iodide) can be detected either by starch (blue) or by PVA (purple-red) in the M-state.22 Using these polyiodide indicators together with pH-indicators revealed separated associated fronts.23,24

■

APPLICATION OF THE DESIGN METHOD IN REAL CHEMICAL SYSTEMS We show step by step how our design method was implemented in experiments and compare the detected spatiotemporal behavior and the types of stationary patterns among the different chemical systems, as well as with the simulation results. We include some practical guidelines and discuss technical and physicochemical factors that one should additionally take into account. Step 1: Selecting pH-Oscillators and Slowly Diffusing Proton Binding Agents

The six pH-systems providing stationary patterns so far are twosubstrate Landolt reactions15 where the first substrate, A−, is sulfite, SO32−. The second substrate C, responsible for proton consumption, is ferrocyanide or alternatively thiosulfate, thiourea, or hydrogen carbonate (Supporting Information Table S1). The reaction paths compiled in reactions R2 and R3 are redox paths, and the oxidant B (either bromate, iodate, or hydrogen peroxide) is actually involved in both paths, except in the case C = hydrogen carbonate, which is a non-redox way to remove protons. Despite the factual sharing of the oxidant, path R3 can be tuned by C independently of path R2 in all these pHsystems. S− was a readily available polyanion, NaPAA, with average 160 carboxylate functions along a single alkyl chain. It is a macromolecule but still small enough to enter the pores of the agarose gel medium. Through the dynamic partitioning between the gel and the CSTR content, the NaPAA concentration in the thin RD medium can be conveniently set by its concentration in the feed solution within a few hours.20,21 3186

DOI: 10.1021/acs.accounts.8b00441 Acc. Chem. Res. 2018, 51, 3183−3190

Article

Accounts of Chemical Research Step 2B: Developing Spatiotemporal Oscillations

When introducing C, the range of the control parameter (here, [H2SO4]0) over which the stability ranges of the F-state and Mstate overlap (spatial bistability) becomes narrower and narrower with [C]0, as seen in simulations. This is usually accompanied by a decrease of the contrast between the two steady states. The onset of spatiotemporal oscillations provides a checkpoint for reaching appropriate negative feedback intensity. The period ranges from 6−11 min (HPSF) to 15−20 min (BSF). The critical value of [C]0 where the two stability limits cross and above which spatiotemporal oscillations develop was similar in the BSF, FIS, and IST systems, while it was much lower in the TuIS and HPSC. The dissimilarity between FIS and TuIS is explained23 by the rate constants of the C + iodine reactions, as proton consumption proceeds via the reduction of the iodine intermediate by C in these systems. Iodine waves (oscillations) were detected in the FIS and IST systems with PVA or starch indicator23,24 but not in the TuIS system, revealing a much lower concentration of the iodine intermediate species. Though calculations with the above skeleton model predict possible bistability between a uniform stationary and the oscillatory state (zone b, Figure 2B), an M/Osc spatial bistability was detected only in the HPSF system so far. This is probably due to the qualitative differences between the boundary conditions in calculations and experiments. However, bistability was evoked to account for a peculiar shrinking dynamics of the oscillating domain as a result of the permanent perturbation induced by the boundary conditions at the rim of the disk OSFR.25 The FIS, TuIS, and BSF systems can develop traveling highpH waves. Spiraling M-state domains (Figure 5a) were observed with the FIS and HPSC reactions.20,21,25 The different combinations of interactions between acid producing (+)fronts and acid consuming (−)fronts resulted in complex dynamics: for example, the repetitive formation and annihilation of an Mstate ring (Figure 5b) in the BSF and BSF−CaEDTA systems.26,27 In the FIS system, repulsion between (−)fronts leads to the formation of more or less long-lived “filaments” of the M-state. Enclosed between filaments, a large oscillatory Mstate domain can trap small F-state “bubbles”, and the dynamics results in the growth and death of F-state spots or in spot splitting (Figure 5c,d).20,21 In the HPSC system, large amplitude extended domains of M-state traveled through quasi-static Mstate filaments, erasing and renewing them at different positions.25 At 30% higher [NaHCO3]0, the domains became smaller and blinked randomly (average period