Collective chemomechanical oscillations in active hydrogels

BZ hydrogel beads loaded with the BZ Ruthenium catalyst [(${\mathrm{Ru}(\mathrm{bpy})}_2(\mathrm{bpyEtAm}){({\mathrm{PF}}_6)}_2$)] are placed into a container of 1 cm in diameter and 5 mm in height, covered with catalyst-free BZ reaction solution, and monitored from above with a digital camera (Fig. 1A). BZ hydrogel beads of radii ranging from $\sim 10$ to 400 $\mathrm{\mu }$m are synthesized following recent procedures (Materials and Methods, A and B) (18, 20). We first create hydrogel beads that are then functionalized with a homemade BZ Ruthenium catalyst. Monodisperse hydrogel beads with tunable concentration of polymer and catalyst are synthesized from a molding technique. Polydispersed hydrogel beads are created from an emulsion polymerization method, allowing fast production of a large number of hydrogel beads with a broad size distribution (21). The resultant radius distribution is fitted with a Gaussian function, with a mean radius of 73 $\mathrm{\mu }$m and a SD of 28 $\mathrm{\mu }$m (SI Appendix, section 2). Here, we use BZ hydrogel beads that have a Ruthenium catalyst concentration $c_{{\mathrm{Ru}}^{2+}}$ ranging from 3 to 30 mM, and a polymer concentration $c_{\mathrm{pol}}$ ranging from 20 to 80 kg m⁻³. The solution of BZ reactant is kept at the following constant concentrations [${\mathrm{H}}_2{\mathrm{SO}}_4$] = 0.4 M, [${\mathrm{NaBrO}}_3$] = 0.08 M, [MA] = 0.06 M, corresponding to a chemical composition widely used for BZ chemomechanical hydrogels. Unless mentioned otherwise, the container is made of a glass slide at the bottom and top, with a PDMS lateral wall of a few millimeters in height and a centimeter in diameter. The solution is kept under pressure with a clamp to avoid evaporation and the creation of carbon dioxide bubbles (Fig. 1A). With this setup, we record the change in color of the catalyst and the radius change of the BZ hydrogel as a function of time (Fig. 1 C and D).

We find that single BZ hydrogel beads do not chemically oscillate when their radius is below a critical radius $R_{\mathrm{c},\exp }=(225\pm 9)\mathrm{\mu }$m (Fig. 1B). For a BZ hydrogel bead with a radius below $R_{\mathrm{c},\exp }$, the catalyst in the BZ hydrogel stays oxidized (Fig. 1C). For a BZ hydrogel bead with a radius above $R_{\mathrm{c},\exp }$, the BZ hydrogel first oscillates with a long period of oscillation and transitions either to a state with a shorter steady period of oscillation or die (Fig. 1D and SI Appendix, section 4). The transition from a long period of oscillation to a shorter and steady period of oscillation occurs through a transient birythmic regime where both frequencies (long and short) are present. We also find that the critical radius $R_{\mathrm{c},\exp }$ does not show a significant variation with the catalyst concentration in the BZ hydrogel bead (SI Appendix, section 3).

To examine how the presence of neighboring BZ hydrogel beads affects their oscillation, we first place a collection of close-packed BZ hydrogel beads with a broad radius distribution in a microfluidic chamber through which we flow the surrounding uncatalyzed BZ solution. We observe that chemical waves of reaction–diffusion propagate in such a divided medium for days (Movie S1). Therefore, BZ hydrogel beads smaller than the critical size $R_{\mathrm{c},\exp }$ can chemically oscillate when surrounded by others. We then assemble quasi-2D layers of BZ hydrogel beads in the chambers we used for the single BZ hydrogel bead experiments (Fig. 1A). Chemical waves of reaction–diffusion are also observed in these 2D media. BZ hydrogel beads as small as $\frac{R_{\mathrm{c},\exp }}{10}\approx 20\mathrm{\mu }$m in radii show sustained chemical oscillations (Fig. 2A). We note that the measured lower size for a BZ hydrogel bead to oscillate is due to the resolution of our optical setup. We also note that quasi-2D layers of BZ hydrogels beads all individually smaller than the critical size for an isolated BZ hydrogel bead to oscillate also support chemical oscillations (Fig. 2B and Movie S2). These experimental results are reminiscent of prior works demonstrating the existence of quorum sensing for BZ microreactors. Oscillating micrometric porous particles loaded with the BZ catalyst dispersed in a solution of uncatalyzed BZ solution have been shown to oscillate in synchrony when strongly stirred (11). The intermediate products generated by the BZ reaction inside the particles diffuse to the outer solution. When stirring is sufficiently strong, the concentration of these mobile chemical species is homogeneous in the uncatalyzed BZ solution. These chemical intermediates mediate the synchronous oscillations of the BZ microreactors. Monolayer made of excitable BZ micrometric particles in a solution of uncatalyzed, unstirred BZ solution has also been found to oscillate as a quorum (12–14). Through numerical simulations, the authors showed that this collective behavior is correlated with a decrease of the activator loss rate averaged over the entire BZ microreactor assembly. In all these previous works and in contrast to our experiments, the BZ microreactors are rigid, hence not mechanically responsive, homogeneous in catalyst concentration and in size, and stay reduced when not oscillating.

To characterize the chemical oscillations in BZ hydrogel beads of various sizes and concentrations of catalyst, we add a few BZ hydrogel beads of specific size and catalyst concentration to a quasi-2D layer of polydispersed BZ hydrogel beads with low catalyst concentration (3 mM) (Fig. 2B). We find that the critical size for a BZ hydrogel bead to oscillate disappears independently of its catalyst concentration (Fig. 2A). Similarly to what we observe with single BZ hydrogel beads, the BZ hydrogel beads experience a transient regime with a long period of oscillation, which then stabilizes to a steady-state regime of oscillation with a constant period (SI Appendix, section 5). The steady-state regime of oscillation, the period of chemical oscillation is found to increase with the catalyst concentration in the BZ hydrogel beads (Fig. 2C).

We now construct a reaction–diffusion theoretical model that qualitatively explains the three most salient features of our experiment: 1) Spherical BZ hydrogel beads smaller than a critical size $R_{\mathrm{c},\exp }$ stay oxidized. Only above a critical radius $R_{\mathrm{c},\exp }$, the BZ hydrogel beads chemically oscillate. 2) Adding neighboring BZ hydrogel beads all smaller than the critical size decreases the critical size for a BZ hydrogel bead to oscillate. 3) In the quorum-sensing regime with heterogeneous BZ hydrogel beads, the steady-state period of oscillation increases with their catalyst concentration.

The physics of an oscillating BZ gel in an uncatalyzed BZ solution is governed by reaction–diffusion processes. The variation over time of the concentration of BZ species $\frac{\partial c_i}{\partial t}$ is due to the diffusion of products in and out of the gel $D\mathrm{\Delta }{c}_i$, with $D$ the molecular diffusion coefficient, and the reaction rate $R(c_i)$ representing the kinetics of the BZ reactions. We simplify this spatially extended model to a point model by assuming the concentration of the BZ species in the BZ hydrogel $c_i$ to be homogeneous and by writing the diffusion term $D\mathrm{\Delta }{c}_i$ as $-{\alpha }_{\mathrm{i}}{c}_i$. Hence,

To justify the mathematical expression of the loss term $-{\alpha }_{\mathrm{i}}{c}_i$, consider the number of moles of a chemical product $i$ inside a BZ gel $N_{i,in}$. This product is diffusing out of the gel with a flux $j_i$ proportional to its concentration gradient, $j_i=-D\mathrm{\nabla }{c}_i\sim D\frac{{\mathrm{dc}}_{\mathrm{i}}}{\mathrm{dx}}\sim {\beta }_i{c}_i$, where ${\beta }_i=\frac{D}{L_i}$, with $L_i$ the typical lengthscale over which the gradient is built across the surface of the bead. The variation of the number of moles of product inside the bead $\mathrm{d}{N}_{i,in}$ due to its loss through its surface in contact with the uncatalyzed BZ solution can be written as $\frac{\mathrm{d}{N}_{i,in}}{\mathrm{d}t}=j_{i}S$, or equivalently as $\frac{\mathrm{d}(V{c}_i)}{\mathrm{d}t}=-{\beta }_i{c}_{i}S$, where we have used that $N_{i,in}=V{c}_i$, with $V$ the volume. Taking the size of a BZ hydrogel bead as constant, the temporal variation of the concentration of product $i$ due to its outward flux becomes $\frac{\mathrm{d}{c}_i}{\mathrm{d}t}=-{\beta }_i\frac{S}{V}{c}_i$. The loss term then acts as an additional effective first-order chemical reaction on the temporal evolution of the chemical concentration. We define ${\alpha }_i={\beta }_i\frac{S}{V}\propto \frac{D}{L_i}\frac{S}{V}$. The lengthscales $L_i$ can be controlled by diffusion, reaction, and reaction-induced convection when the chemicals leaving the hydrogel react in the media outside the hydrogel (18, 22). If we assume that the chemical reactions do not induce convection in the uncatalyzed BZ solution and that these reactions can be modeled as first-order chemical reactions with kinetic constant $K_i$, the lengthscales $L_i$ take the following expression $L_i=\sqrt{\frac{D}{K_i}}$. Hence, ${\alpha }_i$ can be different for each of the mobile species involved in the BZ reaction. In the following, we use the Vanag–Epstein model for the kinetics of the BZ reaction, encoded in the term $R(c_i)$ (23). The Vanag–Epstein (VE) model simplifies the BZ reaction to a model with only four variables, each representing a chemical concentration. The four concentration variables are $c_1$ = [HBrO${}_2$], $c_2$ = [Br${}^{-}$], $c_3$ = [oxidized catalyst] = [Ru${}^{3+}$] and $c_4$ = [Br${}_2$]. HBrO${}_2$, Br${}^{-}$, and Br${}_2$ are chemicals which can diffuse out of the gel, whereas the Ru${}^{3+}$ is bound to the gel. We thus define ${\alpha }_{{\mathrm{Br}}_2}$, ${\alpha }_{{\mathrm{Br}}^{-}}$, and ${\alpha }_{{\mathrm{HBrO}}_2}$ (Materials and Methods, D). We note that while the behavior between different models of the BZ reaction is qualitatively similar, quantitative comparisons between most model simulations of the BZ reaction and experimental results have not been possible (24). Therefore, the quantitative model developed here is intended to reproduce qualitatively the experimental system. For simplicity, we also assume that ${\alpha }_{{\mathrm{Br}}^{-}}={\alpha }_{{\mathrm{Br}}_2}$. Our results do not qualitatively depend on this assumption (SI Appendix, section 7).

With this model, we construct a phase diagram, reported in Fig. 3A, that captures the oscillation state of the BZ hydrogel bead in terms of the rate of loss of the inhibitors ${\alpha }_{{\mathrm{Br}}^{-}}$, ${\alpha }_{{\mathrm{Br}}_2}$ and of the activator ${\alpha }_{{\mathrm{HBrO}}_2}$. We observe 3 regions: At high ${\alpha }_{{\mathrm{HBrO}}_2}$, the BZ hydrogel bead stays reduced, because the activator HBrO${}_2$ diffuses faster than its production rate (red) preventing its autocatalysis to proceed; at high ${\alpha }_{{\mathrm{Br}}^{-}}$ and below certain ${\alpha }_{{\mathrm{HBrO}}_2}$ value, the BZ hydrogel bead stays oxidized, because the inhibitors Br${}^{-}$ and Br${}_2$ diffuse faster than their reaction rates with the ruthenium catalyst (green); otherwise the BZ hydrogel bead oscillates (gray). In the same figure, a change in radius of a BZ hydrogel bead corresponds to a straight line with a slope given by the ratio of the two relevant lengthscales $L_{{\mathrm{HBrO}}_2}$ and $L_{{\mathrm{Br}}^{-}}$. A decrease in the radius of the BZ hydrogel is equivalent to an increase in both ${\alpha }_{{\mathrm{HBrO}}_2}$ and ${\alpha }_{{\mathrm{Br}}^{-}}$. We obtain that the ratio $\frac{L_{{\mathrm{HBrO}}_2}}{L_{{\mathrm{Br}}^{-}}}$ needs to be above 3.6 for a spherical BZ hydrogel to experience a transition from the oxidized state, when small, to the oscillatory state, when its radius is high (Fig. 3 A and C and SI Appendix, section 6). To find a theoretical critical size of $R_{\mathrm{c},\mathrm{theo}}=230\mathrm{\mu }$m consistent with the experimentally observed critical radius $R_{\mathrm{c},\exp }=(225\pm 9)\mathrm{\mu }$m, we use the model with the set of lengthscales L${}_{{\mathrm{HBrO}}_2}$ = 144 $\mathrm{\mu }$m, L${}_{{\mathrm{Br}}^{-}}$ = 34.6 $\mathrm{\mu }$m = L${}_{{\mathrm{Br}}_2}$, resulting in $\frac{L_{{\mathrm{HBrO}}_2}}{L_{{\mathrm{Br}}^{-}}}=4.16$. We verify that the numerical values for L${}_{{\mathrm{HBrO}}_2}$, L${}_{{\mathrm{Br}}^{-}}$, and L${}_{{\mathrm{Br}}_2}$ are in reasonable agreement with the values of these lengthscales calculated from the rate constants of the reactions controlling the concentration profiles of these chemicals outside the BZ hydrogel bead. This is done by numerically integrating the Vanag–Epstein model when the catalyst concentration is set to 0 (Material and Methods, E). This theoretical analysis shows that L${}_{{\mathrm{HBrO}}_2}$, L${}_{{\mathrm{Br}}^{-}}$, and L${}_{{\mathrm{Br}}_2}$ are all smaller than $R_{\mathrm{c}}$. Thus, we find that for a BZ hydrogel bead with a radius $R_c$ at the transition between oscillation and nonoscillation, the reaction rates of the chemicals set the lengthscale over which the chemicals diffuse out of the hydrogel, hence determining the magnitude of their outward flux.

To study the effect of neighboring BZ hydrogel beads on the chemical oscillation, we simulate a line of N hydrogels each affected by the reaction rates of the BZ reaction $R(c_{i,j})$, where $j=1,...,N$ labels the hydrogel bead in the line, and coupled chemically to each other and to the uncatalyzed BZ solution located at the top and side of the line by diffusion. The BZ hydrogel $j$ is coupled to the outside BZ solution by a term $-{\alpha }_i{c}_{i,j}$, an expression we derived earlier for the single BZ hydrogel model. Additionally, the BZ hydrogel $j$ is coupled to its closest neighbors $j-1$ and $j+1$ by the terms $-\beta (c_{i,j}-c_{i,j-1})$ and $-\beta (c_{i,j}-c_{i,j+1})$, respectively. These terms capture the diffusive flux of chemicals between neighboring BZ hydrogel beads (Materials and Methods, F). Solving numerically this system of $4\times \mathrm{N}$ nonlinear equations, we find that the critical size for an assembly of BZ hydrogel beads to oscillate is $R_{\mathrm{c},\mathrm{collective}}=38\mathrm{\mu }$m for all $\beta$ and $N$ values tested. $R_{\mathrm{c},\mathrm{collective}}$ is 6 times smaller than the critical size for the single BZ hydrogel bead to oscillate (Fig. 3C). This result qualitatively reproduces the experimentally observed decrease in the critical size for a BZ hydrogel bead to oscillate in the collective regime. We can recover the same result with a model that considers a 2D layer of BZ hydrogel beads as a BZ hydrogel sheet. The loss terms $-{\alpha }_{\mathrm{i}}$ are proportional to the BZ hydrogel bead surface-to-volume ratio $\frac{S}{V}$. For a sphere of radius R, $\frac{S}{V}=\frac{3}{R}$, whereas for a sheet made of close-packed spherical hydrogels of radius R, $\frac{S}{V}=\frac{1}{2R}$. Hence, ${\frac{S}{V}}_{|\mathrm{sphere}}$ = $6{\frac{S}{V}}_{|\mathrm{sheet}}$. The critical radius for a BZ hydrogel sheet to oscillate is expected to be $R_{\mathrm{c},\mathrm{sheet}}=\frac{1}{6}{R}_{\mathrm{c},\mathrm{sphere}}$ = 38 $\mathrm{\mu }$m, in agreement with the numerical result obtained from the model with N coupled BZ hydrogel beads.

The model with N BZ hydrogel beads also allows simulating the experiments conducted with BZ hydrogel beads containing a heterogeneous catalyst concentration (Fig. 2B). We assume that the parameters ${\alpha }_{i,j}$ are independent of the catalyst concentration. In Fig. 3D, we show that the continuous increase in the oscillation period with catalyst concentration observed experimentally can be qualitatively reproduced in the simulations when the coupling $\beta \le 0.002$s${}^{-1}$. For stronger coupling between neighboring gels, an increase in the oscillation period with the catalyst concentration is also observed, but this increase appears quantized. In this case, $T_{c_{{\mathrm{Ru}}^{2+},high}}=n\times T_{c_{{\mathrm{Ru}}^{2+},low}}$, with $n\in \mathbb{N}$.

To end, we use the collective regime to study the occurrence of chemomechanical oscillations of specific BZ hydrogel beads with radii also below the critical radius $R_{\mathrm{c},\exp }$. These specific BZ hydrogel beads are softer than the rest (bulk modulus $\sim$10 kPa) (18). We find that during an initial transient regime corresponding to a long period of chemical oscillation, such BZ hydrogel beads present mechanical oscillation of a few percent in strain (Movie S3 and SI Appendix, section 8). Their chemical oscillations eventually stabilize to a steady state regime with a faster oscillation period. We have observed that a 150-$\mathrm{\mu }$m-radius BZ hydrogel bead mechanically oscillates in its steady state regime, with an amplitude of around $1\%$ (Fig. 4 and Movie S4). Taken together, these experimental results demonstrate that the collective regime can be leveraged to trigger the chemomechanical oscillation of softer BZ hydrogels.

Using collective effects to change the critical size for a BZ hydrogel bead to chemically oscillate appears as a promising strategy to increase the oscillatory strain. Indeed, as the timescale of swelling of a spherical hydrogel is proportional to the square of its size (25), we hypothesize that decreasing further the size of a BZ hydrogel embedded in this assembly should increase its chemomechanical energy conversion. Future work will consist in testing this idea by tuning the size and environment of BZ hydrogel beads.

We have shown that the collective chemical response of an assembly of soft BZ hydrogel beads triggers their chemomechanical oscillations. A single spherical BZ hydrogel smaller than a critical radius $R_{\mathrm{c},\exp }$ does not chemically oscillate because it stayed oxidized. This is in stark contrast to previous results on the oscillation behavior of catalyst-loaded particles for which the catalyst stays reduced. We then construct a model of the BZ hydrogel bead coupled to its outside environment through diffusion that reproduces this experimental observation. The theory indicates that this size-dependent behavior is due to a competition between the reaction kinetics of the inhibitor of the BZ reaction and its transport to the outer solution. With both experiment and theory, we also show that an assembly of close-packed BZ hydrogel beads with a radius smaller than the previously established critical size for a single BZ hydrogel to oscillate can sustain chemical oscillation. In this collective regime, BZ hydrogel beads, tuned in composition to improve their oscillation strain, are now mechanically active.

These results are important for optimizing the performance of autonomous chemomechanical machines. Indeed, they highlight the role of the size of a BZ hydrogel bead and the role of its environment on its chemomechanical properties. They demonstrate that it is possible to design a multiscale material made with a first lengthscale controlling the mechanical swelling timescale and a second lengthscale dictating the chemical oscillation of the gel. It also shines light on previous experiments made with BZ gels, emulsions or coacervates at the nano/micro scale (26). For example, recent experimental results demonstrate that BZ polymeric shells of around 20 $\mathrm{\mu }$m in diameter undergo buckling unbuckling instabilities (27) and that BZ coacervates or BZ microgels of around 1 $\mathrm{\mu }$m in diameter sustain chemomechanical oscillations (28, 29). Our interpretation of these results is that the oscillations of such small polymeric BZ microreactors might be the result of a collective chemical response. Provided that the density of these BZ microcompartments in the solution is high enough, they can exchange chemicals which can trigger their collective chemomechanical oscillation. Also, in the collective regime, if each BZ hydrogel bead is viewed as a unit cell, our work is a step toward the creation of the recently called “swarmalator,” a model system consisting of a swarm of autonomous chemically communicating chemomechanical oscillators (19, 30–32). Indeed each BZ hydrogel bead consumes chemical energy to communicate with its neighbor and actuate mechanically. The remaining challenge to achieve this goal is to convert its mechanical actuation into motion. The self-propulsion of an autonomous BZ bilayer hydrogel has been experimentally demonstrated by patterning the substrate with a sawtooth design to guide the motion of the BZ bilayer in a specific direction (33). Theoretical work also shows that directed motion can be obtained by adding a gradient of adhesion force between the substrate and the BZ hydrogel (32). In these scenarios, the direction of the BZ gel motion is controlled by the anisotropy of the surface interacting with the BZ gel and is therefore not intrinsic to the BZ gel material itself. As a result, the BZ gel direction cannot adapt to external stimuli such as the chemical gradient generated by neighboring BZ gels. Nonetheless, from these experiments, we can calculate the typical Reynolds number of such BZ gel walkers moving in water. With a speed $u$ = 1 $\mathrm{\mu }$m s${}^{-1}$, a size $L$ = 1 cm, and $\nu ={10}^{-6}$m${}^2$ s${}^{-1}$ the kinematic viscosity of water, we get $Re=\frac{\mathit{uL}}{\nu }\sim 0.01$. At low Reynolds numbers, the self-propulsion is notably hard because breaking time-reversal symmetry is needed to swim (34, 35). The presence of nonlinearity in both mechanics and the BZ chemistry as well as the many degrees of freedom offered by the softness of the hydrogels should facilitate nonreciprocal back and forth motions, hence promoting motility (36–38). Moreover, biology provides examples of organisms that can move at low Reynolds numbers. Those called euglenids harness self-propagating waves directed along their body to move (39, 40). Theoretical work confirms that an active actomyosin network or a BZ chemical reaction on a deformable membrane display autonomous shape changes including a buckling instability that can ultimately lead to motility (41–43). Motile soft robots have also been engineered from 2D membranes powered with a magnetic field (44–46). Therefore, as a first step, 2D shape-changing hydrogels, which have been developed recently, could be made BZ active to generate autonomous folding-unfolding of their structure (47, 48). Further work will then be needed to convert a cyclic shape-changing BZ hydrogel into a self-propelled BZ hydrogel. Taken together, our work opens different directions in the study of the collective behavior of life-like materials, from their locomotion to their self-assembly. Working with a model system such as BZ hydrogels can also help rationalize how multicellular organization emerges as a result of chemical cues and chemomechanical energy conversion at the molecular level.

We thank M. Deforet for critical readings of the manuscript. B.B acknowledges support of the Brandeis Provost research award, the Council for the Arts at Massachusetts Institute of Technology (CAMIT) grant from Massachusetts Institute of Technology (MIT), the financial help from Somerville Art Council, the Maria Zambrano fellowship from the University of Barcelona, and funding from the European Union’s Horizon 2020 research and innovation program under the Marie Sklodowska-Curie grant agreement COFUND Prg Autopo REGION Ile de France (IDF). A.F.-N. thanks MCIN/AEI/10.13039/501100011033/FEDER, UE (grant No. PID2021-122369NB-I00) for financial support.