# Enhancement of chemical oscillations by self-generated convective flows

## Abstract

Chemical feedback loops in fluids can produce not only chemical oscillations, but also density variations that generate solutal buoyancy forces, which in turn initiate fluid flow. Using analytical and computational models, we herein examine how the reaction-induced flows alter the chemical oscillations in a fluid-filled chamber whose top and bottom walls are coated with different enzymes. Due to this chemo-fluidic coupling, the systems form oscillating flow patterns, which combine the characteristic size of the buoyancy-driven convection rolls with the frequency of the chemical oscillations. With changes in the distance between the enzyme-coated walls, the convective flows not only enhance or suppress the chemical oscillations, but also substantially increase the amplitude and frequency of the oscillations and extend the regime of the oscillatory behavior. These design principles can facilitate the development of artificial biochemical networks that act as chemical clocks.

## Introduction

Oscillating chemical reactions in living systems are vital for regulating circadian rhythms, metabolic processes, the transcription of DNA and other crucial biological functions1,2. Within the small-scale confines of a biological cell, molecular diffusion is sufficient to ensure that the reagents are homogeneously mixed and thus, the chemical oscillations depend solely on time and not on space1,2,3,4. On a larger spatial scale, when homogenization due to molecular diffusion cannot be considered instantaneous, the combination of nonlinear chemical reactions and diffusion gives rise to structures that vary in both space and time, such as chemical Turing patterns5 and traveling chemical waves6. This spatio-temporal pattern formation can be properly described by coupled reaction-diffusion equations.

Diffusion is not, however, the only transport mechanism that affects local concentration variations during chemical reactions. If a reaction in solution gives rise to local variations in the fluid density, then the resulting forces (i.e., “solutal buoyancy” forces) will induce the spontaneous motion of the fluid7. A variety of catalytic reactions in solution produce solutal buoyancy forces and thus can generate various complex dynamic behavior8,9,10,11,12. For example, in the presence of the appropriate reactants, surface-bound enzymes can act as “chemical” pumps13, which propel fluids in millimeter-thick chambers. If these enzymes are bound to the inner surfaces of the chamber, the system generates chemicals gradients over a well-defined spatial scale that is determined by the geometry of the chamber. As a result, convective patterns with desirable properties can be efficiently designed. Therefore, these surface-bound enzymatic reactions exhibit advantageous properties for the transport and targeted delivery of submerged cargo14,15,16.

The effects of coupling steady (non-oscillatory) surface reactions and reaction-generated flows was analyzed using linear stability theory17,18. The analysis revealed a similarity in the behavior of the latter system with Rayleigh-Benard convection, where a horizontal liquid layer is heated from below. The applied heat gives rise to density variations in the fluid, resulting in thermal buoyancy19. For fluid convection driven by solutal or thermal buoyancy, stable convective patterns appear when the corresponding buoyancy forces exceed certain critical values. For systems where the fluid flow is coupled to oscillating chemical reactions by the buoyancy or Marangoni mechanisms, a variety of chemical waves and oscillatory fluid flows were observed20,21,22,23,24,25. The latter studies, however, did not explore the rich dynamic behavior and advantages offered by fluidic systems where chemical reactions are coupled to enzymes bound to surfaces, which prove to be beneficial for promoting chemical oscillations and controlling the convective flow.

Here, we model a system where two enzyme-coated surfaces interact through both diffusion and buoyancy-driven convection as reactants are converted to products by the enzymatic reactions. The product of the first enzymatic reactions acts as a promoter for the second reaction. On the other hand, the product of the second reaction acts as an inhibitor for the first. As a result, the system exhibits both positive and negative feedback loops essential for the chemical oscillations. These chemical oscillations are accompanied by oscillating fluid flows that are driven by the solutal buoyancy forces arising in response to variations in the local chemical composition of the solution. The convective flows, in turn, accelerate the transport of chemicals between the enzyme-coated boundaries and, therefore, affect the chemical reactions. We demonstrate that the chemically generated convective flows can, in fact, promote or suppress the chemical oscillations in the system. Moreover, we show that a substantial increase in the amplitude and frequency of the chemical oscillations can occur because of the oscillating convective flows.

## Results

### Theoretical modeling

We consider a solution of A and B chemicals that is confined within an infinitely long horizontal channel of the thickness H; a segment of this channel is shown in Fig. 1. The liquid layer is oriented along the x, y plane, with the vertical z-axis lying perpendicular to the layer. Experimentally, this system could be realized as a “batch reactor” (i.e., no external flow) where the ratio of the height of the chamber to its lateral size is relatively small. As shown in Fig. 1, the bottom wall of the channel is coated with enzyme E1 and the inner, top wall is coated with enzyme E2. In the presence of the substrate S, the enzymes E1 and E2 catalyze the respective transformations of A and B; the corresponding chemical reactions are expressed as $$S\mathop{\longrightarrow}\limits^{{E_1}}A$$ and $$S + A\mathop{\longrightarrow}\limits^{{E_2}}B$$. In addition to these reactions, chemicals A and B undergo deactivation over time; hence, the chemical transformation to some final products is not modeled explicitly. We assume that concentrations of these initial reactants (substrate S) and the final products are time independent6 and the reverse reactions are neglected.

In comparison with chemical reactions promoted by catalysts dissolved in solution1,2,3,4, the spatial separation between the surface-bound enzymes introduces two major advantages. First, binding the enzymes on opposite surfaces (Fig. 1) separates the regions where reactions take place and hence, introduces the time delay required for transporting the participating reactants across the layer. This delay between the production and transformation of the reactants enables chemical oscillations in the system1. The second advantage arises from coupling the surface enzymatic reactions and the fluid motion. The gradients of the reactants introduced by the enzymatic reactions control velocities and the patterns of the generated convective flows.

It is worth noting that the enzymatic reactions shown schematically in Fig. 1 mimic the two-step pathway of the biosynthesis of glutathione26,27, which is a key modulator of the intracellular, chemically reducing environment. Within living cells, the two enzymes are mixed together, and the reaction leads to a steady-state, which varies only due to changes in the environment. The key element of our design is the spatial separation of the two enzymes in a channel. As we show below, the system encompassing the spatially separated enzymes could, under proper conditions, exhibit self-sustaining chemical oscillations, which do not occur within the confines of living cells.

In what follows, we assume that the local density of the solution depends on local concentrations of the reaction products and hence, the solutal buoyancy force due to the spatial variations of local density across the chamber can drive the fluid motion. The solutal buoyancy force is approximately proportional to the gradient of concentration of chemicals in the fluid and is an inherent mechanism for fluid flow. Namely, if one particular reagent is denser than the others and its concentration is locally increased, then locally the fluid becomes denser and sinks downward due to gravity. If the concentration of the latter reagent decreases, the fluid moves upward. These convective fluxes accompany the diffusion of chemicals within the chamber.

For the case shown in Fig. 1, the variations in concentrations CA and CB of the respective reactants A and B change the local density of the solution, Δρ = −ρ0(βACA+  βBCB), where ρ0 is the density of solution at CA = CB= 0. The solutal buoyancy coefficients βj = −(1/ρ0)∂ρ/∂Cj characterize changes in the volume of the solution caused by the dissolved chemicals Cj, j = A, B. With this definition of βj, the less dense products have βj > 0. The corresponding solutal buoyancy force fB = gΔρ acts in the same direction as gravity (g = gez= g(0, 0, 1)) and causes the fluid to move with a velocity u = (ux, uy, uz).

The dynamic behavior of the system is described by a system of equations that includes: the continuity equation

$$\nabla \cdot {\mathbf{u}} = 0,$$
(1)

the Navier–Stokes equation in the Boussinesq approximation28

$$\frac{{\partial {\mathbf{u}}}}{{\partial t}} + \left( {{\mathbf{u}} \cdot \nabla } \right){\mathbf{u}} = - \frac{1}{{\rho _0}}\nabla p + \nu \nabla ^2{\mathbf{u}} + {\mathbf{e}}_zg\left( {\beta ^AC^A + \beta ^BC^B} \right),$$
(2)

where p is pressure, and the reaction-diffusion-advection equations

$$\frac{{\partial C^{{\,}j}}}{{\partial t}} + ({\mathbf{u}} \cdot \nabla )C^j + \gamma ^jC^j = D^j\nabla ^2C^j,\quad j = A,B.$$
(3)

In the above equations, is the spatial gradient operator, ν is the kinematic viscosity, γj is the deactivation (decay) rate constant and Dj is the diffusivity of respective reactants Cj, j = A, B. The expansion coefficients βj couple the chemical and hydrodynamic properties of the system. Namely, if βj= 0 for all j, then there is no buoyancy force and no fluid flow within the channel, i.e., u=0.

The enzymatic reactions occur at the chamber walls and hence, they are introduced through the boundary conditions imposed on Eq. (3). In particular:

$$z = 0:\quad - D^A\frac{{\partial C^A}}{{\partial z}} = k_1\sigma _1F_1(C^B),\\ z = H:\quad - D^A\frac{{\partial C^A}}{{\partial z}} = k_2\sigma _2F_2(C^A),$$
(4)
$$z = 0:\;D^B\frac{{\partial C^B}}{{\partial z}} = 0,\quad z = H:\;D^B\frac{{\partial C^B}}{{\partial z}} = k_2\sigma _2F_2(C^A).$$
(5)

Here, ki is the reaction rate constant and σi is the surface density of enzyme i, where i = 1, 2. The functions F1(CB) and F2(CA) describe the concentration dependence of the inhibited and promoted reactions, respectively, and are chosen to mimic those for the glutathione biosynthesis pathway26,27:

$$F_1(C^B) = \frac{1}{{1 + (C^B/K^B)^{n_1}}},\quad F_2(C^A) = \frac{{(C^A/K^A)^{n_2}}}{{1 + (C^A/K^A)^{n_2}}},$$
(6)

where KB and KA are the respective inhibition and dissociation constants. As seen from Eqs. (5) and (6), the rate of production of chemical A decreases with an increase in the concentration CB (inhibition), whereas an increase in CA increases the rate of production of B until saturation (promotion). Note that the reaction rates in Eqs. (4)–(6) are taken to be dependent on the cooperativity parameters (Hill coefficients) n1,2 > 0. Cooperativity of the enzymatic reactions is known to affect the dynamic regimes that can exist in the system3,4.

Finally, we apply the no-slip boundary conditions for the fluid velocity on the bottom and top surfaces of the chamber:

$${\mathbf{u}}(z = 0) = {\mathbf{u}}(z = H) = {\mathbf{0}}$$
(7)

The problem is analyzed in the dimensionless form described in the Methods. The main dimensionless parameters that govern the problem include the Grashof number, $${\mathrm{Gr}} = \frac{{g\beta ^AH^3C_0}}{{D^{A2}}}$$, which characterizes the magnitude of the buoyancy force, and the dimensionless enzymatic reaction rates $$\mu _\alpha = \frac{{Hk_\alpha \sigma _\alpha }}{{D^AC_0}}$$, α = 1,2, which set the speed of the chemical reactions. Here, C0 is the characteristic concentration of reactants in the system (see Methods). We assume that μ1/μ2 = const and hence, behavior of the system is controlled by the hydrodynamic and chemical processes characterized by the pair (Gr, μ1).

The governing equations, Eqs. (1)–(7), permit the existence of a time independent base state solution with no fluid flow (u=0). The corresponding concentration profiles $${C}_{0}^{A}$$(z) and $${C}_{0}^{B}$$(z) are shown in Fig. 2a. In what follows, we demonstrate that the equilibrium exists only when the parameters Gr and μ1 are below certain critical values, Gr < Grc and μ1 < μc. Otherwise, the equilibrium is unstable, and the instability could result in either the pure chemical oscillations (Gr < Grc, μ1 > μc), convective motion (Gr > Grc, μ1 < μc), or both (Gr > Grc and μ1 > μc). Below, we analyze in detail the onset and interplay of the chemical and hydrodynamic instabilities.

### Oscillatory chemical instability (u = 0)

For sufficiently small expansion coefficients $$(\beta^{{\,}j} < \, \beta_{\mathrm{c}}^{{\,}j})$$, there is no fluid motion (u=0). Oscillations of the concentrations Cj(x, y, z, t) occur when the reaction rate μ1 increases beyond the critical value μc.

The stability curves, which separate the domains of stable and unstable solutions in the plane H − (k1σ1)c, are calculated at the decay constant values of γ = 0.004, 0.001, 0.0005 s−1, and shown in Fig. 2b with solid magenta, dashed green, and dotted azure lines, respectively. The curves corresponding to the different decay constant values are obtained by rescaling the governing parameters as discussed in the Methods. The period of the critical chemical oscillations T = 2π/|ωi| is calculated along the curves in Fig. 2b, and is shown as a function of the layer thickness H in Fig. 2c. As seen in Fig. 2, the linear stability analysis predicts that at a given γ, the chemical instability occurs only within some range of layer thicknesses, Hmin(k1σ1) < H < Hmax(k1σ1), provided that the reaction rate is above a critical value, μ1 > μc. This result demonstrates that the spatial separation between the enzymes is a parameter that controls the chemical oscillations in the system. The results presented in Fig. 2b also imply that larger values of the decay constant γ enable oscillations within narrower channels. In particular, the decay constants of γ = 0.004 s−1 (solid magenta line in Fig. 2b) and γ = 0.0005 s−1 (dotted azure line) enable oscillations in channels with thicknesses around 1 mm and 3–4 mm, respectively. As the oscillation frequency |ωi| = 2π/T is inversely proportional to the period of oscillation T, the decrease in T observed for increasing values of γ (shown with azure, green, and magenta lines in Fig. 2c) also means an increase in the frequency of the chemical oscillations.

Having identified the most favorable parameters (H, μc(H)) for the growth of chemical oscillations at q = 0, we now fix the thickness H, and examine the oscillatory instability with respect to the non-zero wave number q = 2π/λ. The stability curve (k1σ1)c(q) plotted in Fig. 3a shows that the most rapidly growing oscillations have zero wave number, q* = 0, thereby indicating that the chemical instability, which exists at u=0, is of the long-wavelength type.

### Buoyancy-driven monotonic instability

At the low reaction rates, μ1 < μc, when the chemical oscillations do not exist, the stability of the state of mechanical and chemical equilibrium can be broken by increasing the chemical expansion coefficients beyond their critical values $$\beta^{{\,}j}\, > \, \beta_{\mathrm{c}}^{{\,}j}$$. The onset of steady convection is determined by solving the eigenvalue problem, Eqs. (17)–(20), at ωr = 0. Specifically, we found only the monotonic instability for which ωr = 0, ωi = 0, and a critical Raleigh number, Rac, is determined at a given wave number of perturbation, q. (The solutal Raleigh number is defined as $${\mathrm{Ra}} = \frac{{g\beta ^AH^3C_0}}{{\nu D^A}}$$, see Methods).

Figure 3b shows the numerically obtained stability boundary, which separates the domain of steady convection (u0) above the curve from the steady-state with u=0 below the curve. The buoyancy-driven instability is plotted in the dimensionless variables as is commonly done in the literature28. As seen in Fig. 3b, convection exists at Ra > Rac ≈ 1447, which is first achieved at the wave number of q* = 2π/λ* = 3. At Ra > Rac, the steady convective motion is possible within a window of wavenumbers around q* (Fig. 3b), and the perturbations with q = q* have the largest growth rate. It is worth noting that the critical Rayleigh number Rac depends on the reaction rates μα as discussed in the Methods.

Thus, the analysis of linear stability predicts the existence of instabilities having two different mechanisms. The first mechanism is the chemical long-wavelength oscillatory instability at μ1 >  μc(q) with the boundary μc(q) defined by the imaginary growth rate ωr = 0 and ωi ≠ 0. The latter instability exists at subcritical Rayleigh numbers, Ra < Rac, and the minimum value μc(q*) is reached for zero wave number q* = 0 (Fig. 3a). The instability arises as a result of the interaction between the diffusion and nonlinear chemical reactions of the chemicals A and B.

The second instability is a monotonic buoyancy-driven short-wavelength instability at Ra > Rac(q) with the boundary Rac(q) defined by the real growth rate ωr = 0 and ωi = 0. The instability arises in response to the variations in the fluid density due to the local chemical composition of the solution. The denser fluid volumes sink to the bottom while the less dense rise to the top of the channel. This instability exists for the subcritical reaction rate μ1 < μc, and the boundary reaches its minimum Rac(q*) at q* = 3 as seen in Fig. 3b.

The results of linear stability analysis are summarized in the diagram in Fig. 3c. Values of the critical reaction rate μc(q*), indicating transitions from the base state to the chemical oscillations, are shown with a vertical dashed green line while the transitions to the convective motion (u0) at Rac(q*) are marked by a horizontal solid line. The nonlinear regimes beyond the stability domain with the boundaries μ1 = μc and Ra = Rac, are discussed in the next section.

### Two-dimensional (2D) regimes with subcritical Rayleigh numbers

To investigate the system far from the stability boundaries, we simulate the system’s behavior through the numerical solution14 of Eqs. (9)–(13) in a 2D cell, where 0 ≤ x ≤ L, 0 ≤ z ≤ H, and the periodic boundary conditions are applied in the x direction. The cell length L is set equal to the length of the critical perturbation that gives rise to convection, L =  λ*, which is determined through linear theory. Simulations in a longer domain L = 5λ* are discussed in the Supplementary Note 2. As initial conditions, we use the uniform spatial distribution of reactants Cj(x, z, t = 0) = rj, where 0 ≤ rj ≤ 1 is a random number, and take u(x, z, 0) =  0 (no fluid motion). The parameters μ1 and Ra characterize the external “stimuli”, which can drive the chemical oscillations or fluid motion or both.

We first focus on the region encompassing subcritical values of the Rayleigh number (Ra < Rac) where u=0. When μ1 < μc in this region, the stimuli are weak, and the linear stability analysis indicates that the system is in a stable steady-state (Fig. 3c) with no fluid motion and a time independent distribution of the concentration profiles $$C_0^{{\,}j}(z)$$ across the layer (Fig. 2a). On the other hand, at the supercritical reaction rate μ1 > μc (Fig. 3c), the linear stability analysis predicts that the concentrations of chemicals A and B, Cj(z,t), exhibit temporal oscillations that are spatially uniform along the layer, i.e., independent of x.

The simulation results on the onset of chemical oscillations are in agreement with the predictions from the linear stability analysis (Figs. 2 and 3a, c). Namely, the plots in Fig. 4 show that the chemical oscillations (with no convective flow u=0) exist at the supercritical reaction rates μ1 >  μ1c and subcritical Rayleigh numbers Ra < Rac. Figure 4a displays the temporal variations in the concentrations CA (red) and CB (blue) that are averaged over the z coordinate, 〈C j(z, t)〉z. Figure 4b shows the spatial profiles of CA (red) and CB (blue) along the z direction at the moments of time corresponding to the maximal (solid lines) and minimal (dashed lines) values exhibited by the averaged concentrations in the course of oscillations (see Fig. 4a).

Figure 4c reveals that the amplitudes of the chemical oscillations, $${A}_{\mathrm{chem}}^{j}$$ = A, B, as functions of the reaction rate μ1 exhibit the behavior characteristic for a supercritical bifurcation. The amplitudes were calculated as $${A}_{\mathrm{chem}}^{j}$$ = max〈C j(z,t)〉z − min〈C j(z,t)〉z and plotted with solid black squares connected by red (for A) and blue (for B) solid lines. Notably, $${A}_{\mathrm{chem}}^{j}$$ =  0 at μ1 ≤ μc, with μc being marked by the vertical dashed green line. Above the bifurcation point, μc, the oscillation amplitudes exhibit an increase, which is approximately proportional to the square root of the distance from the bifurcation point, $${A}_{\mathrm{chem}}^{j}\propto$$ (μ1 − μc)1/2. Figure 4c also shows that the values of average concentrations smoothly continue through the critical point μc from the sub- into the supercritical region, as can be seen by the open circles connected by red (for CA) and blue (for CB) lines. The average values 〈C j(z,t)〉z,t were obtained by averaging over the z coordinate and over a period of time much longer than the period of oscillation. Finally, Fig. 4d shows that the period of oscillations, T, increases with an increase in the reaction rate k1σ1.

### Two-dimensional regimes with supercritical Rayleigh numbers

In the supercritical region, Ra > Rac, the solutal buoyancy effects can give rise to convective flow. It is reasonable to expect that at the subcritical reaction rate μ1 < μc (see the stability curve in Fig. 3b), the convection should begin as a result of the instability of the base state, whereas at μ1 > μc, both the convection and chemical oscillations should be present (the stability of the uniform oscillations is not studied in the above calculations). Below, we discuss the details of the system’s behavior at Ra > Rac as obtained by the numerical simulations.

For sufficiently small reaction rates μ1 < μc and supercritical Ra > Rac, the fluid motion emerges in the form of steady convective rolls as predicted by the linear stability theory (Fig. 3b, c). The flow field in the established regime is shown with black arrows in Fig. 5a. The amplitude of the fluid motion is characterized by the square of the maximal fluid velocity, $$u_{\mathrm{max}}^2(t)$$, which is plotted (red dots in Fig. 5b) as a function of the expansion coefficient βA (~Ra). On a computational mesh that contains Nx and Nz nodes in the x and z directions, respectively, we calculate the maximal fluid velocity as $$u_{\mathrm{max} }^2(t)={\mathrm{max}}\{u_{x}^2(i,j,t)+u_z^2(i,j,t), 1\le i \le N_x,1\le j \le N_z\}$$. To characterize the oscillatory regimes, we plot the value $$\overline {u_{\mathrm{max} }^2} = \frac{1}{2}\left( {{\mathrm{max} } \{ u_{\mathrm{max} }^2(t)\} + \min \{ u_{\mathrm{max} }^2(t)\} } \right)$$, which is a measure of the maximal fluid velocity during one cycle of a stable oscillation. Near the bifurcation point, the amplitude of the convective rolls increases with an increase in the expansion coefficient as approximately $$u_{\mathrm{max}}^2(t)\propto \beta^A-\beta_{\mathrm{c}}^A\, (u_{\mathrm{max}}^2(t) \propto {\mathrm{Ra}}-{\mathrm{Ra}}_{\mathrm{c}})$$. The latter behavior is presented by the magenta line (k1σ1 = 0.76 μmol m−2 s−1) in Fig. 5b. The amplitude of the fluid motion for the steady convective rolls is indicated by the red dots in Fig. 5b.

For reaction rates close enough to the critical value, μ1 ~ μc (shown with a dashed red line in Fig. 5b), and at Ra > Rac(μ1), the steady convective rolls become unstable, and the new stable regime combines the concentration oscillations with oscillations of the convective rolls. The amplitude of the fluid motion is presented by the black squares in Fig. 5b.

The oscillations in the convective flow have a spatial period of L = 2, which is imposed by the boundary conditions and is consistent with critical wavelength λ* of the buoyancy-driven instability. The oscillation frequency of the oscillatory convective flow is found to be Ω ≈ ωi, i.e., close to the frequency of the chemical oscillations. (Oscillations in a longer domain L = 5λ* are presented in Supplementary Fig. 1 and Supplementary Movie 2). The oscillatory convective flow regime exists within a window of the solutal expansion coefficients $$\beta_{\mathrm{c}}^A\, < \, \beta^A \, < \, \beta_2^A$$ (Rac < Ra < Ra2) as seen in Fig. 5b. An increase in the expansion coefficient βA (~Ra) leads to an increase in the amplitude of the oscillating flow, $$u_{\mathrm{max}}^2(t)$$, and Fig. 5b shows that $$u_{\mathrm{max}}^2(t)$$ increases much faster than the linearly growing amplitude for the rolls $$u_{\mathrm{max}}^2(t) \propto \beta^A-\beta_{\mathrm{c}}^A$$. This increase in the amplitude of the convective oscillations occurs approximately at the onset of the instability for the steady convective rolls ($$\beta \sim \beta_{\mathrm{c}}^{A}$$), reaches the maximal values, and gradually becomes suppressed by the rolls at the transition point $$\beta_{2}^{A}$$. The red dashed line in Fig. 5b indicates the critical value μc, which separates the amplitudes of the convective oscillations for the subcritical μ1 < μc and supercritical μ1 > μc reaction rates.

It is worth noting that in the presence of convective flows, the results of the 2D simulations reveal that convective chemical oscillations exist in both the subcritical μ1 < μc (below the dashed red line in Fig. 5b) and supercritical μ1 > μc reaction rates. In other words, the presence of the convective flow enlarges the domain of chemical oscillations, which do not exist in the subcritical region μ1 < μc in the absence of the fluid flow (see Fig. 4c). The existence of chemical oscillations at μ1 < μc in the presence of convective flows cannot be captured by the linear stability analysis (see Fig. 3c). This nonlinear effect can only be detected and studied by numerical solution of the system of nonlinear equations, Eqs. (9)–(11).

The temporal and spatial behaviors of the oscillatory convection regime are shown in Fig. 6. Oscillations of the concentrations 〈CAx,z and 〈CBx,z averaged across the area of the domain are shown in Fig. 6a with the red and blue lines, respectively. While the chemical oscillations at u=0 exhibit a harmonic-like behavior (Fig. 4a), the concentration oscillations are seen to be strongly nonlinear in the presence of the convective flow (Fig. 6a). The difference between the mean values of the concentrations 〈CAx,z,t and the concentrations 〈CBx,z,t in Fig. 6a is smaller than that for purely chemical oscillations in Fig. 4a. Oscillations of the maximal value of the fluid velocity in the middle of the domain are shown in Fig. 6b. The peaks in the values of maximal velocity coincide in time with the corresponding maxima of the reactant concentration 〈CAx,z, which is coupled to the fluid motion through the expansion coefficient βA. The structure of the fluid flow oscillates between the pattern in Fig. 6c, which corresponds to the peak values of the fluid velocity, and the almost motionless state in Fig. 6d, which corresponds to the minimal velocity values. Flow structures for the entire oscillation period are described in Supplementary Note 3 and presented in Supplementary Fig. 2 and Supplementary Movie 1.

Notably, the continuation of the amplitude $$u_{\max }^2(t) \propto \beta^A-\beta_{\mathrm{c}}^A$$ of the convective rolls to the instability onset $$\beta^{A}=\beta_{\mathrm{c}}^{A}$$ suggests that the oscillatory regime originates approximately at the same critical value $$\beta^{A}=\beta_{\mathrm{c}}^{A}$$. We hypothesize that the oscillatory regime may develop as a result of a secondary instability of the steady convective rolls. Some limited evidence supporting this hypothesis comes from the results of simulation at k1σ1 = 7.6 and 3.8 μmol m−2 s−1 (cyan and blue lines in Fig. 5a). Namely, the results indicate that for parameters βA ≈ 0.01 M−1 close to the instability onset $$\beta_{\mathrm{c}}^{A}$$, the convective rolls (red dots) occur for smaller values of βA than the convective oscillations (black squares). As the equilibration of the steady regimes near the bifurcation point is very slow, further investigations are needed to clarify the origin of the observed behavior.

The effect of the convective flow on the chemical oscillations is evident upon comparing the amplitude and frequency of the concentration oscillations during the oscillatory convection with those at no flow (u=0) shown in Fig. 4. As faster fluid flows enable faster transport of the reactants between the enzyme-coated walls of the chamber, both the amplitude and frequency of the convective oscillations increase with an increase in the expansion coefficient βA (~Ra) as demonstrated in Fig. 7. Figure 7a, b show the amplitude, $${A}_{\mathrm{chem}}^{j}$$ = max〈CAx,z − min〈CAx,z, and period of oscillation, T, of the concentration CA as functions of the expansion coefficient βA within the range $$\beta_{\mathrm{c}}^{A} < \, \beta ^A < \, \beta_{2}^{A}$$ where the oscillations exist. As seen in Fig. 7, the amplitude, $${A}_{\mathrm{chem}}^{j}$$, increases and the period of oscillation, T, decreases with an increase in the reaction rate k1σ1. For comparison, the horizontal green dashed line shows the oscillation amplitude with no flow (u=0) taken from Fig. 4c for k1σ1 = 11.4 μmol m−2s−1. Also, at u=0, the period of chemical oscillations is an increasing function of k1σ1 as seen in Fig. 4d. Notably, the value of k1σ1 affects the range of $$\beta^{A}, \beta_{\mathrm{c}}^{A} \, (\mu_1) \, < \, \beta^{A} < \, \beta_2^A \, (\mu_1)$$, where the oscillation exists. The observed increase in the amplitude and frequency of oscillations is due to the intensification of the reactant transport across the layer by the reaction-generated convective flows.

## Discussion

We developed and analyzed a theoretical model of a system where the chemical oscillations and density-driven convective flows are enabled by the enzymes bound to the top and bottom walls of the horizontal channel filled with the solution of reactants. In contrast to typical models of nonlinear chemical dynamics29,30, our model introduces the nonlinear reactions through the boundary conditions to the reaction-diffusion equations. This study extends our understanding of how chemical and hydrodynamic degrees of freedom could interact within biochemical and synthetic systems.

The distance between the enzyme-coated walls was shown to control both the chemical oscillations and the speed of the generated fluid flows. We found that the oscillating convective flows combined the characteristic size of the buoyancy-driven convection rolls with the frequency of the chemical oscillations.

It is worth noting that the interaction between the short-wavelength Turing patterns and long-wavelength chemical oscillations, which leads to formation oscillatory patterns, occurs within the Brusselator model31. In the Brusselator model, the short and long-wavelength modes have the same origin rooted in diffusion and nonlinear chemical kinetics. In contrast, the regime of the convective oscillations described within our model arises as a result of the interaction between short-wavelength convective rolls, which have a hydrodynamic origin, and long-wavelength chemical oscillations governed by diffusion and chemical kinetics.

We showed that in the millimeter-thick fluidic channels, the solutal buoyancy-driven convection can substantially increase the amplitude and frequency of the oscillations and enlarge the domain of their existence. The latter behavior is due to the acceleration of transport of the reactants between the enzyme-covered walls by the convective fluxes. Sufficiently fast convective flows, however, were found to suppress the chemical oscillations because the oscillations do not exist in a well-mixed solution.

The minimal model developed here involves only two intermediate reactants and a buoyancy force that arises due to the presence one of these reactants. Notably, the solutal buoyancy force can arise for a broad variety of chemical reactions where molecules of the reactants and reaction products occupy different volume in the solution. The magnitude of the reported effects is determined by the velocities of the chemically generated fluid flows, which advect reactants and directly modify chemical fluxes. Both the magnitude of the buoyancy force and the geometry of the container determine the velocity of the generated flows. For example, in sufficiently narrow channels, the fluid motion is suppressed by the no-slip boundary conditions, which ensure zero velocities on the walls. Qualitatively, the flow speed is characterized by the magnitude of Grashof number Gr, which multiplies the buoyancy force term in Eq. (10) and is proportional to both the cube of the channel thickness, H3, and the characteristic density variation, ρ0βAC0. Thus, velocities of the generated convective flows are more sensitive to the chamber thickness than to the magnitude of the solutal buoyancy force. Therefore, for given reactions and given buoyancy forces, the thickness of the chamber, H, can be adjusted to generate fluid flows with velocities that maximize the reported effects of the flow on chemical oscillations.

The design principles based on the mechanisms described here could facilitate the development of artificial biochemical networks that involve chemical clocks. Note that the period of oscillations in biochemical networks1,2, where the components are uniformly mixed in the solution, is typically on the order of hours. A design that employs surface-bound enzymes can accelerate the oscillations through the generated convective fluid flows.

Finally, we note that there are other mechanisms besides the solutal buoyancy force that could be used to generate fluid motion. In case of exothermic reactions, the generated heat could give rise to an additional buoyancy force $$\propto g\beta_{T}\Delta T$$ proportional to the expansion coefficient βT = −(1/ρ0)∂ρ/∂T and the local temperature increase ΔT. We ignored this effect because the typical values of solutal expansion coefficients (βC = 0.01–0.1 M−1) are larger than the values of the thermal buoyancy coefficients (βT ~ 10−4 K−1)14,15,16,32. The oscillatory chemical reactions could also be coupled to the fluid motion through the thermo- or soluto-capillary22,23,24 (in the case of free liquid surfaces) and diffusiophoretic33 effects specified by the corresponding boundary conditions. In all these cases, we expect the generated flows to modify the chemical oscillations, either suppressing or amplifying these oscillations.

## Methods

### Dimensionless variables describing the system

To simplify the equations, we introduce dimensionless parameters. The width of the channel H is chosen as the characteristic length scale. The two time scales, H2/ν, which characterizes the viscosity of the solution, and H2/DA, which is related to the reactant diffusivity, are related by the Schmidt number $${\mathrm{Sc}} = \frac{\nu }{{D^A}}$$ ~ 103 (Supplementary Table 1). We define the diffusion time scale in order to introduce the dimensionless variables describing the system

$${\mathbf{r}}^{\prime} = \frac{{\mathbf{r}}}{H},\quad t^{\prime} = \frac{{D^At}}{{H^2}},\quad {\mathbf{u}}^{\prime} = \frac{{H{\mathbf{u}}}}{{D^A}},\quad p^{\prime} = \frac{{H^2}}{{(D^A)^2\rho _0}}p,\quad C^{\prime j} = \frac{{C^j}}{{C_0}},\\ K^{\prime j} = \frac{{K^j}}{{C_0}},\quad \mu^{\prime}_\alpha = \frac{{Hk_\alpha \sigma _\alpha }}{{D^AC_0}},$$
(8)

where j = A, B, α = 1, 2 and C0 is the characteristic concentration of reactants. We set C0 = 2.1 M, which is the maximal value of the base state concentration calculated for parameters in Supplementary Table 1, i.e., $$C_0\equiv C_0^{A}(z=0)$$. After dropping the prime symbols, the equations governing behavior of the system, Eqs. (1)–(3), are reduced to the following dimensionless form:

$$\nabla \cdot {\mathbf{u}} = 0$$
(9)
$$\frac{{\partial {\mathbf{u}}}}{{\partial t}} + \left( {{\mathbf{u}} \cdot \nabla } \right){\mathbf{u}} = - \nabla p + {\mathrm{Sc}}\nabla ^2{\mathbf{u}} + {\mathbf{e}}_z{\mathrm{Gr}}\left( {C^A + \delta C^B} \right)$$
(10)
$$\frac{{\partial C^j}}{{\partial t}} + ({\mathbf{u}} \cdot \nabla )C^j + \Gamma ^jC^j = \xi ^j\nabla ^2C^j,\quad j = A,B,$$
(11)

where Γj and ξj are dimensionless parameters defined further below. The corresponding boundary conditions, Eqs. (4), (5), (7), can be written as

$$z = 0: - \frac{{\partial C^A}}{{\partial z}} = \mu _1F_1(C^B),\quad z = 1: - \frac{{\partial C^A}}{{\partial z}} = \mu _2F_2(C^A), \\ z = 0:\frac{{\partial C^B}}{{\partial z}} = 0,\quad z = 1:\frac{{\partial C^B}}{{\partial z}} = \frac{{\mu _2}}{{\xi ^B}}F_2(C^A),$$
(12)
$${\mathbf{u}}(z = 0) = {\mathbf{u}}(z = 1) = {\mathbf{0}}.$$
(13)

As seen in Eqs. (10)–(12), behavior of the system is controlled by the following dimensionless parameters:

$${\mathrm{Gr}} = \frac{{g\beta ^AH^3C_0}}{{D^{A2}}},\quad \delta = \frac{{\beta ^B}}{{\beta ^A}},\quad \Gamma ^j = \frac{{H^2\gamma ^j}}{{D^A}},\quad {\mathrm{Sc}} = \frac{\nu }{{D^A}},\\ \xi ^j = \left\{ \begin{array}{l}1,\quad j = A\\ \frac{{D^B}}{{D^A}},j = B,\end{array} \right.\quad {\mathrm{Ra}} = \frac{{{\mathrm{Gr}}}}{{{\mathrm{Sc}}}}$$
(14)

The dimensionless parameters include the solutal Grashof number, Gr, which describes the ratio between the buoyancy and viscous forces, the Schmidt number, Sc, which characterizes the ratio between the diffusive and viscous time scales, and the solutal Rayleigh number, Ra, which is the ratio between the solutal Grashof and Schmidt numbers. Note that the functions Fα (α = 1, 2) now depend on the dimensionless variables defined by Eq. (8).

The dimensionless formulation of the problem implies that any solution to the problem, Eqs. (9)–(13), for a given set of parameters H, γ, kασα, and βj is also a solution for the rescaled parameters $$H^*= \zeta H, \gamma^*=\gamma \zeta^{-2}, k_{\alpha}^*\sigma_{\alpha}^*=k_{\alpha}\sigma_{\alpha}\zeta ^{-1}$$, and βj* = βjζ−3, where ζ is the dimensionless scaling parameter. This commonality in solutions holds because the rescaling preserves the dimensionless parameters Γj, μα, and Gr in Eqs. (9)–(13) that describe the system’s behavior.

To simplify the analysis, we reduce the number of independent dimensionless parameters by assuming that DA = DB= D and γA = γB= γ, so that Γj= Γ = H2γ/D and ξj = 1 for j = A,B. We assume that the fluid motion arises only due to variations in the concentration CA, i.e., that βA ≠ 0, βB≠ 0 and hence, δ = 0. We also set the Hill coefficients to n1 = n2 = 3. Finally, we take KA = KB = K, and the dimensionless reaction rates to be μ1/μ2 = const ≈ 0.0288 with μ1 being an independent variable. The values of the physical properties of the chemical solution are given in Supplementary Table 1.

### Base state solution

The horizontally homogeneous base state of the system is characterized by the mechanical equilibrium, at which the concentrations are stationary, ∂Cj/∂t = 0, and there is no fluid motion, u = 0. The base state is controlled by the expansion coefficient βA and described by the equations:

$${\frac{{\partial p_0}}{{\partial z}} = {\mathrm{Gr}}\,C_0^A(z)} \\ {d^2C_0^j/dz^2 = \Gamma \,C_0^j}.$$
(15)

The concentrations $$C_0^j$$, j = A, B, depend only on the coordinate across the layer, z, and can be presented as

$$C_0^j(z) = a_0^j{\mathrm{e}}^{\Lambda z} + b_0^j{\mathrm{e}}^{ - \Lambda z},$$
(16)

where Λ = Γ1/2. In the above equation, the coefficients $$a_0^j$$ and $$b_0^j$$ depend on the reaction rates μα, α = 1, 2, and are found from the boundary conditions given by Eq. (12). The profiles $$C_0^A(z)$$ and $$C_0^B(z)$$ calculated for H = 2 mm and k1σ1 = 9.85 μmol m−2 s−1 are shown by the respective red and blue lines in Fig. 2a.

The hydrodynamic and chemical processes are governed by the Grashof number Gr and enzymatic reaction rates μα. The system behavior is characterized by the pair (Gr, μ1) since we assumed that μ1/μ2 = const.

Note that the Grashof number, Gr, and the dimensionless reaction rate μ1 are defined through a certain characteristic concentration C0 (see the Eqs. (14) and (8)). As usual for problems involving convection, the characteristic concentration is defined as the difference in the chemical concentration $$C_0^A(z)$$ across the layer, $$C_0 = C_0^A(1)-C_0^A(0)$$, calculated in the stationary state, Eq. (16), at βA = 0, i.e., when convection does not exist. As the value of C0 depends on the reaction rate μ1, the Grashof number is an implicit function of the reaction rate, Gr(C0(μ1)).

### Linear stability analysis

We study the stability of the base state, Eqs. (15), (16), by introducing small perturbations u = 0 + u1(x, y, z, t), p = p0(z) + p1(x, y, z, t) and $$C^{{\,}j}=C_0^{{\,}j}(z)+ C_1^{{\,}j}(x,y,z,t)$$ with normal modes $$\{ {u_{z1},p_1,C_1^j} \} = \left\{ {f\left( z \right),p\left( z \right),c^j\left( z \right)} \right\}{\mathbf{e}}^{(\omega t + iq_xx + iq_yy)}$$, which are periodic along the layer. The dynamics of the perturbations is described by the following boundary value problem (Supplementary Note 1):

$$- \omega \left( {\frac{{\partial ^2f}}{{\partial z^2}} - q^2f} \right) + {\mathrm{Sc}}\left( {\frac{{\partial ^4f}}{{\partial z^4}} - 2q^2\frac{{\partial ^2f}}{{\partial z^2}} + q^4f} \right) = q^2{\mathrm{Gr}}\,c^A$$
(17)
$$\frac{{\partial ^2c^j}}{{\partial z^2}} - \left( {\omega + \Gamma ^j + q^2} \right)c^j = f\frac{{\partial C_0^j}}{{\partial z}},\quad j = A,B$$
(18)
$$f(z = 0) = \left. {\frac{{\partial f}}{{\partial z}}} \right|_{z = 0} = 0,\quad f(z = 1) = \left. {\frac{{\partial f}}{{\partial z}}} \right|_{z = 1} = 0$$
(19)
$$z = 0: - \frac{{\partial c^A}}{{\partial z}} = \mu _1F^{\prime}_1(C_0^B)c^B,\quad z = 1: - \frac{{\partial c^A}}{{\partial z}} = \mu _2F^{\prime}_2(C_0^A)c^A, \\ z = 0: - \frac{{\partial c^B}}{{\partial z}} = 0,\quad z = 1:\frac{{\partial c^B}}{{\partial z}} = \mu _2F^{\prime}_2(C_0^A)c^A,$$
(20)

Here, the wave number $$q = \frac{{2\pi }}{\lambda } = \sqrt {q_x^2 + q_y^2}$$ and complex frequency ω = ωr + i characterize the size λ and time evolution of the perturbation. In addition, here $$i = \sqrt { - 1}$$, and the prime in $$F_{\alpha}^{\prime}$$, α =  1, 2, denotes the derivative with respect to Cj. Eqs. (17)–(20) are solved numerically using the shooting method34.

### Oscillatory chemical instability

For sufficiently small expansion coefficients $$(\beta^j \, < \, \beta_{\mathrm{c}}^j)$$, there is no fluid motion (u=0) and the problem (Eqs. (1720)) reduces to the equation d2cj/dz2 = Ωj2cj, which describes oscillatory chemical instability and has the solution

$$c^j(z) = a^j{\mathrm{e}}^{\Omega z} + b^j{\mathrm{e}}^{ - \Omega z},\quad \Omega ^2 = \Gamma + \omega + q^2,\quad j = A,B$$
(21)

subject to the boundary conditions in Eq. (20).

The instability boundary, which separates the parameters enabling chemical oscillations from those supporting the stationary distribution of chemicals A and B, is obtained by solving Eqs. (20) and (21) for the imaginary eigenvalue ω = i at q = 0. The Hill coefficients n1,2 in Eq. (6) determine the response of the enzyme to concentration variations. Consistent with other studies where oscillatory behavior is described by the Hill function (as in Eq. (6))3,4, we find no oscillatory behavior for n1,2 < 3 regardless of the values of other parameters.

### Buoyancy-driven monotonic instability

The onset of steady convection is determined by solving the eigenvalue problem, Eqs. (17)–(20), at ωr = 0. Specifically, we found only the monotonic instability for which ωr = 0, ωi = 0, and a critical Raleigh number, Rac, is determined at a given wave number of perturbation, q.

## Data availability

The data that support the findings of this study are available in the main text and supplementary information. Additional information is available from the corresponding author upon request.

## Code availability

The code is available from the corresponding author upon reasonable request.

## References

1. 1.

Novak, B. & Tyson, J. Design principles of biochemical oscillators. Mol. Cell Biol. 9, 981 (2008).

2. 2.

Lim, W. A., Lee, C. M. & Tang, C. Design principles of regulatory networks: searching for the molecular algorithms of the cell. Mol. Cell 49, 202 (2013).

3. 3.

Elowitz, M. B. & Leibler, S. A synthetic oscillatory network of transcriptional regulators. Nature 403, 335 (2000).

4. 4.

Shum, H., Yashin, V. V. & Balazs, A. C. Self-assembly of microcapsules regulated via the repressilator signaling network. Soft Matter 11, 3542 (2015).

5. 5.

Turing, A. M. The chemical basis of morphogenesis. Philos. Trans. R. Soc. Lond., Ser. B 237, 37 (1952).

6. 6.

Prigogine, I. & Lefever, R. Symmetry breaking instabilities in dissipative systems. II. J. Chem. Phys. 48, 1695 (1968).

7. 7.

Pojman, J. A. & Epstein, I. R. Convective effects on chemical waves. 1. Mechanisms and stability criteria. J. Phys. Chem. 94, 4966 (1990).

8. 8.

Rongy, L., Trevelyan, P. M. J. & De Wit, A. Dynamics of A + B → C reaction front in the presence of buoyancy-driven convection. Phys. Rev. Lett. 101, 084503 (2008).

9. 9.

Almarcha, C., Trevelyan, P. M. J., Grosfils, P. & De Wit, A. Chemically driven hydrodynamic instabilities. Phys. Rev. Lett. 104, 044501 (2010).

10. 10.

Sebestikova, L. & Hauser, M. J. B. Buoyancy-driven convection may switch between reactive states in three-dimensional chemical waves. Phys. Rev. E 85, 036303 (2012).

11. 11.

Rossi, F., Budroni, M. A., Marchettini, N. & Carballido-Landeira, J. Segmented waves in a reaction-diffusion-convection system. Chaos 22, 037109 (2012).

12. 12.

Zhang, Y., Tsitkov, S. & Hess, H. Complex dynamics in a two-enzyme reaction network with substrate competition. Nat. Catal. 1, 276 (2018).

13. 13.

Sengupta, S. et al. Self-powered enzyme micropumps. Nat. Chem. 6, 415 (2014).

14. 14.

Shklyaev, O. E., Shum, H., Sen, A. & Balazs, A. C. Harnessing surface-bound enzymatic reactions to organize microcapsules in solution. Sci. Adv. 2, e1501835 (2016).

15. 15.

Das, S. et al. Harnessing catalytic pumps for directional delivery of microparticles in microchambers. Nat. Commun. 8, 14384 (2017).

16. 16.

Maiti, S., Shklyaev, O. E., Balazs, A. C. & Sen, A. Self-organization of fluids in a multienzymatic pump system. Langmuir 35, 3724 (2019).

17. 17.

Bdzil, J. & Frisch, H. L. Chemical instabilities. II. Chemical surface reactions and hydrodynamic instability. Phys. Fluids 14, 475 (1971).

18. 18.

Mc Taggart, C. L. & Straughan, B. Chemical surface reactions and nonlinear stability by the method of energy. SIAM J. Math. Anal. 17, 342 (1986).

19. 19.

Rayleigh, L. On convection currents in a horizontal layer of fluid, when the higher temperature is on the under side. Philos. Mag. 32, 529 (1916).

20. 20.

Matthiessen, K. & Muller, S. C. Global flow waves in chemically induced convection. Phys. Rev. E 52, 492 (1995).

21. 21.

Kitahata, H., Aihara, R., Magome, N. & Yoshikawa, K. Convective and periodic motion driven by chemical wave. J. Chem. Phys. 116, 5666 (2002).

22. 22.

Diewald, M., Matthiessen, K., Muller, S. C. & Brand, H. R. Oscillatory hydrodynamic flow due to concentration dependence of surface tension. Phys. Rev. Lett. 77, 4466 (1996).

23. 23.

Miike, H., Miura, K., Nomura, A. & Sakurai, T. Flow waves of hierarchical pattern formation induced by chemical waves: The birth, growth and death of hydrodynamic structures. Phys. D. 239, 808 (2010).

24. 24.

De Wit, A., Eckert, K. & Kalliadasis, S. Introduction to the focus issue: chemo-hydrodynamic patterns and instabilities. Chaos 22, 037101 (2012).

25. 25.

Budroni, M. A. et al. Control of chemical chaos through medium viscosity in a batch ferroin-catalyzed Belousov-Zhabotinsky reaction. Phys. Chem. Chem. Phys. 19, 32235 (2017).

26. 26.

Jez, J. M., Cahoon, R. E. & Chen, S. Arabidopsis thaliana glutamate-cysteine ligase. J. Biol. Chem. 279, 33463 (2004).

27. 27.

Jez, J. M. & Cahoon, R. E. Kinetic mechanism of glutathione synthetase from Arabidopsis thaliana. J. Biol. Chem. 279, 42726 (2004).

28. 28.

Chandrasekhar, S. Hydrodynamic and Hydro Magnetic Stability (Oxford, Clarendon, 1961).

29. 29.

Scott, S. K. Oscillations, Waves, and Chaos in Chemical Kinetics (Oxford, New York, 1994).

30. 30.

Epstein, I. R. & Pojman, J. A. An Introduction to Nonlinear Chemical Dynamics: Oscillations, Waves, Patterns, and Chaos (Oxford, New York, 1998).

31. 31.

Yang, L., Zhabotinsky, A. M. & Epstein, I. R. Stable squares and other oscillatory Turing patterns in a reaction-diffusion model. Phys. Rev. Lett. 92, 198303 (2004).

32. 32.

Valdez, L., Shum, H., Ortiz-Rivera, I., Balazs, A. C. & Sen, A. Solutal and thermal buoyancy effects in self-powered phosphatase micropumps. Soft Matter 13, 2800 (2017).

33. 33.

Derjaguin, B. V., Sidorenkov, G. P., Zubashchenko, E. A. & Kiseleva, E. V. Kinetic phenomena in boundary films of liquids. Kolloidn. Zh. 9, 335 (1947).

34. 34.

Stoer, J. & Burlisch, R. Introduction to Numerical Analysis (New York, Springer-Verlag, 1980).

## Acknowledgements

This work was supported as part of the Center for Bio-Inspired Energy Science, an Energy Frontier Research Center funded by the US Department of Energy, Office of Science, Basic Energy Sciences under Award DE-SC0000989.

## Author information

Authors

### Contributions

O.E.S. developed the theory and performed the analysis; V.V.Y. identified parameters crucial for the effect; S.I.S. suggested the reaction scheme; A.C.B. organized the work and analyzed the data. All authors contributed in manuscript preparation and data interpretation discussions.

### Corresponding author

Correspondence to Anna C. Balazs.

## Ethics declarations

### Competing interests

The authors declare no competing interests.

Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Supplementary information

### 42005_2020_341_MOESM3_ESM.avi

Supplementary Movie 1

### 42005_2020_341_MOESM4_ESM.avi

Supplementary Movie 2

## Rights and permissions

Reprints and Permissions

Shklyaev, O.E., Yashin, V., Stupp, S.I. et al. Enhancement of chemical oscillations by self-generated convective flows. Commun Phys 3, 70 (2020). https://doi.org/10.1038/s42005-020-0341-3
{"extraction_info": {"found_math": true, "script_math_tex": 0, "script_math_asciimath": 0, "math_annotations": 0, "math_alttext": 0, "mathml": 0, "mathjax_tag": 0, "mathjax_inline_tex": 0, "mathjax_display_tex": 2, "mathjax_asciimath": 0, "img_math": 0, "codecogs_latex": 0, "wp_latex": 0, "mimetex.cgi": 0, "/images/math/codecogs": 0, "mathtex.cgi": 0, "katex": 0, "math-container": 0, "wp-katex-eq": 0, "align": 0, "equation": 0, "x-ck12": 0, "texerror": 0, "math_score": 0.9416126608848572, "perplexity": 1088.7473726558612}, "config": {"markdown_headings": true, "markdown_code": true, "boilerplate_config": {"ratio_threshold": 0.18, "absolute_threshold": 10, "end_threshold": 15, "enable": true}, "remove_buttons": true, "remove_image_figures": true, "remove_link_clusters": true, "table_config": {"min_rows": 2, "min_cols": 3, "format": "plain"}, "remove_chinese": true, "remove_edit_buttons": true, "extract_latex": true}, "warc_path": "s3://commoncrawl/crawl-data/CC-MAIN-2020-24/segments/1590348511950.89/warc/CC-MAIN-20200606062649-20200606092649-00363.warc.gz"}