## Abstract

We developed an efficient model for responsive gels that captures large-scale, two-dimensional (2D) deformations and chemical reactions within a swollen polymer network. The 2D calculations allowed us to probe not only volume changes but also changes in sample shape. By focusing on gels undergoing the oscillatory Belousov-Zhabotinsky reaction, we observed traveling waves of local swelling that form a rich variety of dynamic patterns and give rise to distinctive oscillations in the gel's shape. The observed patterns depend critically on the gel's dimensions. The approach provides a useful computational tool for probing the dynamics of chemomechanical processes and uncovering morphological transformations in responsive gels.

For a synthetic material to perform sustained mechanical work, it must undergo large-scale, periodic changes in volume or shape. Polymeric gels constitute optimal candidates for use as soft active materials, because controlled modulations of the surrounding solvent can lead to significant, rhythmic expansion and contraction of the gel (*1*, *2*). Consequently, such oscillating gels could be used as microactuators (*3*) for pulsatile drug delivery (*4*). In addition to their practical utility, oscillating gels provide an ideal medium for investigating nonlinear dynamical phenomena that can arise from a coupling of mechanical and chemical energy. For example, researchers have isolated scenarios where the initial swelling and deswelling of a chemoresponsive gel exerts feedback on a nonoscillatory chemical reaction and thereby drives the entire system into a regular, oscillatory regime (*4*–*8*). Theoretical models for oscillating gels (*2*, *7*–*9*) have yielded substantial insight into mechanisms that can produce periodic pulsations. Prior calculations (*2*, *7*–*9*), however, were carried out in one dimension (e.g., the sample was assumed to be spherically symmetric), and thus only volumetric changes of the sample could be probed. To capture shape changes, models must describe the gel deformation in at least two dimensions. By encompassing additional degrees of freedom, 2D (or 3D) models can enhance our fundamental understanding of the interplay between the finite deformations of a responsive medium and the nonlinear chemical dynamics and open up the possibility of uncovering new morphological transitions (*7*).

We developed an approach for simulating chemoresponsive gels that can exhibit not only large variations in volume but also alterations in shape. Through this approach, we modeled oscillating gels undergoing the Belousov-Zhabotinsky (BZ) reaction (*10*–*12*) and showed that the ensuing 2D pattern formation depends on the dimensions of the sample. BZ gels are unique because the polymer network can expand and contract periodically without external stimuli (*11*, *12*). This autonomous, self-oscillatory behavior is due to a ruthenium catalyst, which is covalently bonded to the polymers (*11*, *12*). The BZ reaction generates a periodic oxidation and reduction of the anchored metal ion, and these chemical oscillations induce the rhythmic swelling-deswelling in the gel. By probing morphological changes in BZ gels, we can establish design criteria for creating autonomous smallscale devices, which perform sustained work until the reagents in the host solution are consumed and can be simply refueled by replenishing these solutes.

Our approach was inspired by the lattice spring model (LSM), where a material is represented by a network of interconnected Hookean springs (*13*); however, we made substantial modifications to create a gel lattice spring model, or gLSM, which captures the gel's large-scale deformations, interaction with the surrounding solvent, and response to chemical reactions. We started with a 2D lattice, which corresponds to a thin film or a 3D sample where the height of the material remains constant. In the initial undeformed system, this layer height is *H*_{0}, and the length of each element in the square lattice is Δ (fig. S1A). The spatial position of every node in this lattice is uniquely characterized by the vector **m** = *k* e_{1} + *l* e_{2} or, equivalently, by thepairofintegers(*k,l*). The elements are labeled by the pair **m** = (*k*, *l*) corresponding to the node in the lower left corner of the element. Upon deformation, the vectors **m** merely label the nodes and the elements, whereas the vectors *x*_{m} give the actual positions of the nodes in the laboratory coordinate system (fig. S1A). The total energy of the deformed system, *E*, can be expressed as a function of the coordinates of the nodes, or *E* = *U*({*x*_{m}}), where *U*({*x*_{m})} = *H*_{0}Δ^{2}Σ_{m}*u*_{m} is the sum of the energies of the deformed rectangular elements.

To capture the nonlinear elasticity of the gels, we must express the energy density, *u*, in terms of the invariants of the strain tensor, *I*_{i} (*14*). For purely 2D deformation, the invariants of interest are (*14*, *15*) and where *d*_{i} = ∂** x**/∂

*X*

_{i}for

*i*values of 1 and 2, and λ

_{⊥}and

**e**

_{3}are the degree of swelling and the unit vector in the direction perpendicular to the deformation, respectively. The energy density

*u*=

*u*(

*I*

_{1},

*I*

_{3}) is a sum of two distinct contributions (

*8*). The first,

*u*

_{el}, describes the rubber elasticity of the cross-linked polymer chains and can be expressed as (

*16*) (1) where

*T*is temperature (measured in energy units) and

*c*

_{0}is the cross-link density, that is, the number density of elastic strands in the undeformed polymer network. In the undeformed state,

*u*

_{el}equals 0 because

*I*

_{1}= 3 and

*I*

_{3}= 1. The second contribution to the energy density,

*u*

_{FH}, describes the interaction between the polymer and solvent units: (2) where

*v*

_{0}is the volume of a monomer within the chains,

*f*

_{FH}is the Flory-Huggins free energy density (

*17*), and

*f*

_{FH}(ϕ) = (1 – ϕ) ln(1 – ϕ) + χ

_{FH}ϕ(1 – ϕ), with χ

_{FH}being the polymer-solvent interaction parameter. The local volume fraction of polymer in the gel, ϕ, depends on the volume fraction of polymer in the undeformed

_{1}state, ϕ

_{0}, and the local deformation as .

The force acting on a node **m** can now be obtained as *F*_{m} = – ∂*U*/∂*x*_{m}, which yields (3) Here, *P*_{m} is the pressure in element **m**, with *P*_{m} = π_{m} + *c*_{m}*T*/2. In the latter equation, π_{m} is the Flory-Huggins osmotic pressure (*17*) in element **m**: (4) where is the volume fraction of polymer in element **m**. Similarly, *c*_{m} is the local value of the cross-link density and is equal to . Because of the high friction between the polymer and solvent, inertial effects can be neglected, and the node velocity is proportional to the force acting on that site (*18*) so that (5) The mobility is inversely to the polymer-solvent friction coefficient, ζ, where for a good solvent (*18*).

Equation 3 allows us to interpret the internal gel forces as a sum of two contributions: (i) the spring-like forces acting between nodes and (ii) the forces that originate from the pressure between the adjacent elements. A graphical representation of the equation for the force acting on node **m** =(*k, l*) is shown in fig. S1B. Note that one can readily introduce physical and chemical heterogeneities at different locations in the sample and thereby analyze how the local variations affect global properties.

This gLSM can be readily generalized to incorporate the BZ reaction. In addition to the local volume fraction of polymer, ϕ_{m}, and the cross-link density, *c*_{m}, each element is now characterized by the dimensionless concentrations of the following (*9*, *19*): reactant dissolved in the solvent, *u*_{m}, and oxidized metal-ion catalyst attached to the chains, *v*_{m}. The gel dynamics are coupled to the BZ reaction via the pressure, *P*_{m}, as (6) where χ* is the coupling constant. When χ* > 0, the gel can respond to the chemical oscillations of the BZ reaction by undergoing spatial oscillations; when χ* = 0, the gel is nonresponsive to the chemical oscillations. For each element of the gLSM, the equations of chemical kinetics are formulated in terms of the mole fractions of the reagents, and Because the metal-ion catalyst is chemically linked to the polymer, the local chemical balance equation includes only the redox reaction: (7) where ϵ is a constant. The other reactant undergoes both reaction and diffusion: (8) The first term on the right-hand side of Eq. 8 describes the diffusive transport between the neighboring elements. **J**_{i}(**m**) is the flux of the dissolved reagent across edge *i* of the element **m**, with ξ_{i}(**m**) being the vector associated with this edge (fig. S1A). Both the transport of this reagent due to the polymer-solvent interdiffusion and the self-diffusion of the reagent in the solvent contribute to the flux J_{i}(**m**). *G* and *F* characterize the BZ reaction in the gel through the modified (*9*) Oregonator model (*19*): (9) (10) where the constants *f*, ϵ, and *q* are the model parameters; a parameter accounting for the light sensitivity of the BZ reaction can be incorporated into Eq. 9 as outlined in (*20*). Equations 5, 7, and 8 are solved numerically to obtain ϕ_{m}, *u*_{m}, and *v*_{m} as a function of time. We used the VODE code for solving large systems of stiff ordinary differential equations (*21*). The gel is freely swelling (i.e., no external forces act on the boundary nodes), and *u* equals 0 in the outer solution. Because the model is purely 2D, there are no transport processes in the third direction. The reaction rate functions *F* and *G* depend on ϕ, which depends on the local deformations; thus, the mechanical response of the system could affect the chemical kinetics of the system, as will be seen below.

The density plots in Fig. 1 show the temporal evolution of ϕ and the concentration of *v* for a particular portion of the gel. Specifically, we focused on the marked cross sections of a square sample (Fig. 1A) and plotted ϕ and *v* along these lines as a function of time. To highlight the difference between the responsive and nonresponsive gels, we first consider a system in which χ* = 0; that is, the polymer dynamics are not coupled to the oscillatory chemical reaction. The appearance of repeating, straight stripes in the *v* density plot (Fig. 1B) indicates that this variable exhibits periodic temporal oscillations (as dictated by the chosen values of *f*, ϵ, and *q*); however, the uniform coloring of the ϕ density plot reveals that the oscillatory changes in the catalyst have no effect on the chains.

When χ* > 0, we found dramatically different behavior. The coupling of the BZ reaction to the elasticity of the gel gives rise to traveling waves within the polymeric network, as seen by the waves of yellow in ϕ (Fig. 1C). The fact that these bent curves point to the left indicates that the waves originate from the center and propagate to the edges of the sample as time increases. In addition to this wave of local swelling, the BZ reaction modifies the shape of the responsive sample, as can be seen by the periodic variations in the sample's width. In this square example, similar oscillations appear in the vertical direction, causing the gel to undergo a symmetric, rhythmic expansion and contraction. The mechanical oscillations of the gel actually modify the BZ reaction, as can be seen by comparing *v* in Fig. 1, B and C. [It is worth noting that, for BZ reactions in solution, convective motion can lead to a feedback that alters the chemical reaction (*22*).]

We found that the dimensions of the sample have a notable effect on the dynamic patterns and shape changes that occur within the gel. Figure 2 shows snapshots in time for three examples: The sample width, *w*, is fixed at 40 lattice sites, whereas the height, *h*, is set to 10, 20, and 30 sites. The outer edges are marked in black to highlight the morphological transformations. For *h* = 10, the chemical waves lead to a swelling at the center of the slab; the swelling then propagates to the edges of the sample to yield a dog bone–shaped object. The entire process repeats, with the structure returning to the shape in column I, row A (IA). In this case, the waves emanate from the center and spread to the outer edges. This behavior also characterizes the patterns for the *h* = 30 sample, shown in column III. A swelling rapidly propagates to the edges, leaving an unswollen core in IIIF. The system then repeats the cycle, restarting with the configuration in IIIA.

The observed patterns are, however, much more intriguing when *h* = 20 (*23*). The waves appear to originate from two foci (Fig. 2, column II), merge in the center, and then decorate the edges, so that system resembles two back-to-back Cs. The outer edges of the Cs join to form the pattern in IIF, and the oscillations continue, returning to the pattern in IIA.

To further illustrate the distinct nature of the patterns for *h* = 20, we plotted the structural evolution of this system in both the lateral (Fig. 3) and vertical (Fig. 4) directions. Along the a-b cross section (Fig. 3A), the spatiotemporal pattern appears as repeating chevrons that point to the right (Fig. 3B). This pattern is indicative of waves that originate from the edges and propagate toward the center of the sample as time moves forward. In Fig. 3C, which shows the behavior along the c-d cross section, each bright line consists of two chevrons, corresponding to two waves that propagate away from the points of origin. When these waves reach the sample's center, they form a wave that propagates toward the edges in the vertical direction (Fig. 4B). If the vertical cross section passes through one of the two foci (Fig. 2 II), then the evolution exhibits a standing-wave pattern (Fig. 4C). Lastly, near the outer edges of the sample, the wave propagates inward again (Fig. 4D). Thus, the dynamic behavior for *h* =20 is much more complex than those seen for *h* =10 and *h* = 30, where the waves simply progress toward the edges, generating patterns similar to those in Fig. 1C. Although the aspect ratio *h*/*w* plays a role, further studies are needed to quantify how this value affects pattern formation.

To establish the relevant time and length scales, we noted that the unit of time in the dimensionless Eqs. 6 and 7 is *T*_{0}, the shortest time scale of the BZ reaction in solution (*19*, *24*). On the basis of rate constants in solution (*25*) and the typical concentrations in gels (*26*), we estimated *T*_{0} ∼ 1 s. The period of time in Fig. 2 is δ*t* = 100 *T*_{0} or on the order of 1 to 2 min. The unit of length is *L*_{0} = (*D _{u}T_{0}*)

^{½}, where

*D*

_{u}= 2 × 10

^{–5}cm

^{2}s

^{–1}is the diffusion constant of the reactant

*u*in solution (

*27*), so the estimated length scale is Δ =

*L*

_{0}∼ 40 μm. These time and length scales are experimentally accessible, allowing the predictions to be readily tested.

By applying our gLSM technique to responsive gels undergoing the BZ reaction, we uncovered a rich variety of dynamic patterns and distinctive shape changes that depend on the dimensions of the sample. Additionally, within this system, chemical and mechanical behavior are highly coupled and mutually responsive. Overall, the studies further facilitate the design of devices that harness chemomechanical energy conversion to exhibit self-sustained rhythmic action in multiple directions.