Enhancement of chemical oscillations by self-generated convective flows

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. Oscillating chemical reactions are ubiquitous in living systems. Here, the authors propose a mathematical model to study how chemically-driven density gradients produced by enzymatic reactions can trigger hydrodynamic instabilities coupling with chemical oscillations.

O scillating chemical reactions in living systems are vital for regulating circadian rhythms, metabolic processes, the transcription of DNA and other crucial biological functions 1,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 space [1][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 patterns 5 and traveling chemical waves 6 . 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 fluid 7 . A variety of catalytic reactions in solution produce solutal buoyancy forces and thus can generate various complex dynamic behavior [8][9][10][11][12] . For example, in the presence of the appropriate reactants, surfacebound enzymes can act as "chemical" pumps 13 , 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 cargo [14][15][16] .
The effects of coupling steady (non-oscillatory) surface reactions and reaction-generated flows was analyzed using linear stability theory 17,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 buoyancy 19 . 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 observed [20][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 E 1 and the inner, top wall is coated with enzyme E 2 . In the presence of the substrate S, the enzymes E 1 and E 2 catalyze the respective transformations of A and B; the corresponding chemical reactions are expressed as S À!
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 independent 6 and the reverse reactions are neglected.
In comparison with chemical reactions promoted by catalysts dissolved in solution [1][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 system 1 . 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 glutathione 26,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 Fig. 1 Schematic of the system. Segment of an infinite horizontal liquid layer, of thickness H, confined between two enzyme-coated plates (red and blue lines at z = 0 and z = H) promoting chemical transformations of reactants A (red) and B (blue). (x, z) is the coordinate frame and g is the gravity vector.
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 C A and C B of the respective reactants A and B change the local density of the solution, Δρ = −ρ 0 (β A C A + β B C B ), where ρ 0 is the density of solution at C A = C B = 0. The solutal buoyancy coefficients β j = −(1/ρ 0 )∂ρ/∂C j characterize changes in the volume of the solution caused by the dissolved chemicals C j , j = A, B. With this definition of β j , the less dense products have β j > 0. The corresponding solutal buoyancy force f B = gΔρ acts in the same direction as gravity (g = ge z = g(0, 0, 1)) and causes the fluid to move with a velocity u = (u x , u y , u z ).
The dynamic behavior of the system is described by a system of equations that includes: the continuity equation the Navier-Stokes equation in the Boussinesq approximation 28 where p is pressure, and the reaction-diffusion-advection equations In the above equations, ∇ is the spatial gradient operator, ν is the kinematic viscosity, γ j is the deactivation (decay) rate constant and D j is the diffusivity of respective reactants C j , 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: Here, k i is the reaction rate constant and σ i is the surface density of enzyme i, where i = 1, 2. The functions F 1 (C B ) and F 2 (C A ) describe the concentration dependence of the inhibited and promoted reactions, respectively, and are chosen to mimic those for the glutathione biosynthesis pathway 26,27 : where K B and K A 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 C B (inhibition), whereas an increase in C A 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) n 1,2 > 0. Cooperativity of the enzymatic reactions is known to affect the dynamic regimes that can exist in the system 3,4 .
Finally, we apply the no-slip boundary conditions for the fluid velocity on the bottom and top surfaces of the chamber: The problem is analyzed in the dimensionless form described in the Methods. The main dimensionless parameters that govern the problem include the Grashof number, Gr ¼ gβ A H 3 C 0 D A2 , which characterizes the magnitude of the buoyancy force, and the dimensionless enzymatic reaction rates μ α ¼ Hk α σ α D A C 0 , α = 1,2, which set the speed of the chemical reactions. Here, C 0 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 A 0 (z) and C B 0 (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 < Gr c and μ 1 < μ c . Otherwise, the equilibrium is unstable, and the instability could result in either the pure chemical oscillations (Gr < Gr c , μ 1 > μ c ), convective motion (Gr > Gr c , μ 1 < μ c ), or both (Gr > Gr c 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 ðβ j < β j c Þ, there is no fluid motion (u = 0). Oscillations of the concentrations C j (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 − (k 1 σ 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, H min (k 1 σ 1 ) < H < H max (k 1 σ 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 (k 1 σ 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 longwavelength 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 β j > β j c . 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, Ra c , is determined at a given wave number of perturbation, q. (The solutal Raleigh number is defined as Ra ¼ gβ A H 3 C 0 νD A , see Methods). Figure 3b shows the numerically obtained stability boundary, which separates the domain of steady convection (u ≠ 0) 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 literature 28 . As seen in Fig. 3b, convection exists at Ra > Ra c ≈ 1447, which is first achieved at the wave number of q* = 2π/λ* = 3. At Ra > Ra c , 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 Ra c 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 < Ra c , 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 shortwavelength instability at Ra > Ra c (q) with the boundary Ra c (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 Ra c (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 (u ≠ 0) at Ra c (q*) are marked by a horizontal solid line. The nonlinear regimes beyond the stability domain with the boundaries μ 1 = μ c and Ra = Ra c , 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 solution 14 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 Fig. 3 Stability boundaries as functions of the wave number. a Critical reaction rate as a function of the wave number q produces a boundary of a chemical long-wavelength oscillatory instability, which separate regions of the base state C 0 j (z) (below the dashed curve) and chemical oscillations C 0 j (x,z,t) (above) in the absence of the fluid flow, i.e., u = 0. The minimal critical reaction rate (k 1 σ 1 ) c (q) is achieved at the zero wave number q* = 0 (red dot). The layer thickness is H = 2 mm. b The instability curve Ra c (q) formed by critical Rayleigh number as a function of the wave number q separates regions of the base state (u = 0, C 0 j (z)) from the convective fluid flows (u ≠ 0). The minimal critical value Ra c (q c ) = 1447 is achieved at q* = 3 (red square). c Schematic phase-diagram of the regimes obtained from the linear stability analysis as a function of the critical reaction rate μ c (q* = 0) and Rayleigh number Ra c (q* = 3). initial conditions, we use the uniform spatial distribution of reactants C j (x, z, t = 0) = r j , where 0 ≤ r j ≤ 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 < Ra c ) 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 j 0 ð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, C j (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 < Ra c . Figure 4a displays the temporal variations in the concentrations C A (red) and C B (blue) that are averaged over the z coordinate, 〈C j (z, t)〉 z . Figure 4b shows the spatial profiles of C A (red) and C B (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 j chem = A, B, as functions of the reaction rate μ 1 exhibit the behavior characteristic for a supercritical bifurcation. The amplitudes were calculated as A j chem = 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 j chem = 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 j chem / (μ 1 − μ c ) 1/2 . Figure 4c also shows that the values of average concentrations smoothly continue through the critical point μ c from the subinto the supercritical region, as can be seen by the open circles connected by red (for C A ) and blue (for C B ) 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 k 1 σ 1 .

Two-dimensional regimes with supercritical Rayleigh numbers.
In the supercritical region, Ra > Ra c , 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 > Ra c as obtained by the numerical simulations.
For sufficiently small reaction rates μ 1 < μ c and supercritical Ra > Ra c , 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 2 max ðtÞ, which is plotted (red dots in Fig. 5b) as a function of the expansion coefficient β A (~Ra). On a computational mesh that contains N x and N z nodes in the x and z directions, respectively, we calculate the maximal fluid velocity as u 2 max ðtÞ ¼ maxfu 2 x ði; j; tÞ þ u 2 z ði; j; tÞ; 1 ≤ i ≤ N x ; 1 ≤ j ≤ N z g. To characterize the oscillatory regimes, we plot the value u 2 max ¼ 1 2 maxfu 2 max ðtÞg þ minfu 2 max ðtÞg À Á , 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 2 max ðtÞ / β A À β A c ðu 2 max ðtÞ / Ra À Ra c Þ. The latter behavior is presented by the magenta line (k 1 σ 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 > Ra c (μ 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 β A c < β A < β A 2 (Ra c < Ra < Ra 2 ) as seen in Fig. 5b. An increase in the expansion  Fig. 5b shows that u 2 max ðtÞ increases much faster than the linearly growing amplitude for the rolls u 2 max ðtÞ / β A À β A c . This increase in the amplitude of the convective oscillations occurs approximately at the onset of the instability for the steady convective rolls (β $ β A c ), reaches the maximal values, and gradually becomes suppressed by the rolls at the transition point β A 2 . 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 〈C A 〉 x, z and 〈C B 〉 x,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 〈C A 〉 x,z,t and the concentrations 〈C B 〉 x,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 〈C A 〉 x,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 2 max ðtÞ / β A À β A c of the convective rolls to the instability onset β A ¼ β A c suggests that the oscillatory regime originates approximately at the same critical value β A ¼ β A c . 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 k 1 σ 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 β A c , 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 enzymecoated 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 j chem = max〈C A 〉 x,z − min 〈C A 〉 x,z , and period of oscillation, T, of the concentration C A as functions of the expansion coefficient β A within the range β A c < β A < β A 2 where the oscillations exist. As seen in Fig. 7, the amplitude, A j chem , increases and the period of oscillation, T, decreases with an increase in the reaction rate k 1 σ 1 . For comparison, the horizontal green dashed line shows the oscillation amplitude with no flow (u = 0) taken from Fig. 4c for k 1 σ 1 = 11.4 μmol m −2 s −1 . Also, at u = 0, the period of chemical oscillations is an increasing function of k 1 σ 1 as seen in Fig. 4d. Notably, the value of k 1 σ 1 affects the range of 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 dynamics 29,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 shortwavelength Turing patterns and long-wavelength chemical oscillations, which leads to formation oscillatory patterns, occurs within the Brusselator model 31 . 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 enzymecovered 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, H 3 , and the characteristic density variation, ρ 0 β A C 0 . 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 networks 1,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 / gβ T Δ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~1 0 −4 K −1 ) [14][15][16]32 . The oscillatory chemical reactions could also be coupled to the fluid motion through the thermo-or soluto-capillary [22][23][24] (in the case of free liquid surfaces) and diffusiophoretic 33 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, H 2 /ν, which characterizes the viscosity of the solution, and H 2 /D A , which is related to the reactant diffusivity, are related by the Schmidt number Sc ¼ ν D A~1 0 3 (Supplementary Table 1). We define the diffusion time scale in order to introduce the dimensionless variables describing the system where j = A, B, α = 1, 2 and C 0 is the characteristic concentration of reactants. We set C 0 = 2.1 M, which is the maximal value of the base state concentration calculated for parameters in Supplementary Table 1, i.e., C 0 C A 0 ðz ¼ 0Þ. After dropping the prime symbols, the equations governing behavior of the system, Eqs. (1)- (3), are reduced to the following dimensionless form: where Γ j and ξ j are dimensionless parameters defined further below. The corresponding boundary conditions, Eqs. (4), (5), (7), can be written as As seen in Eqs. (10)- (12), behavior of the system is controlled by the following dimensionless parameters: 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).
To simplify the analysis, we reduce the number of independent dimensionless parameters by assuming that D A = D B = D and γ A = γ B = γ, so that Γ j = Γ = H 2 γ/D and ξ j = 1 for j = A,B. We assume that the fluid motion arises only due to variations in the concentration C A , i.e., that β A ≠ 0, β B ≠ 0 and hence, δ = 0. We also set the Hill coefficients to n 1 = n 2 = 3. Finally, we take K A = K B = 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, ∂C j /∂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: The concentrations C j 0 , j = A, B, depend only on the coordinate across the layer, z, and can be presented as where Λ = Γ 1/2 . In the above equation, the coefficients a j 0 and b j 0 depend on the reaction rates μ α , α = 1, 2, and are found from the boundary conditions given by Eq. (12). The profiles C A 0 ðzÞ and C B 0 ðzÞ calculated for H = 2 mm and k 1 σ 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 C 0 (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 A 0 ðzÞ across the layer, C 0 ¼ C A 0 ð1Þ À C A 0 ð0Þ, calculated in the stationary state, Eq. (16), at β A = 0, i.e., when convection does not exist. As the value of C 0 depends on the reaction rate μ 1 , the Grashof number is an implicit function of the reaction rate, Gr (C 0 (μ 1 )).
Linear stability analysis. We study the stability of the base state, Eqs. (15), (16), by introducing small perturbations u = 0 + u 1 (x, y, z, t), p = p 0 (z) + p 1 (x, y, z, t) and C j ¼ C j 0 ðzÞ þ C j 1 ðx; y; z; tÞ with normal modes fu z1 ; p 1 ; C j 1 g ¼ f z ð Þ; p z ð Þ; f c j z ð Þge ðωtþiq x xþiq y yÞ , which are periodic along the layer. The dynamics of the perturbations is described by the following boundary value problem (Supplementary Note 1): Àω Here, the wave number q ¼ 2π λ ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi q 2 x þ q 2 y q and complex frequency ω = ω r + iω i characterize the size λ and time evolution of the perturbation. In addition, here i ¼ ffiffiffiffiffiffi À1 p , and the prime in F 0 α , α = 1, 2, denotes the derivative with respect to C j . Eqs. (17)- (20) are solved numerically using the shooting method 34 .
Oscillatory chemical instability. For sufficiently small expansion coefficients ðβ j < β j c Þ, there is no fluid motion (u = 0) and the problem (Eqs. (17)(18)(19)(20)) reduces to the equation d 2 c j /dz 2 = Ω j 2 c j , which describes oscillatory chemical instability and has the solution c j ðzÞ ¼ a j e Ωz þ b j e ÀΩz ; 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ω i at q = 0. The Hill coefficients n 1,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 n 1,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, Ra c , 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.