Modeling autonomously oscillating chemo-responsive gels

Progress in Polymer Science 35 (2010) 155–173

Contents lists available at ScienceDirect

Progress in Polymer Science journal homepage: www.elsevier.com/locate/ppolysci

Modeling autonomously oscillating chemo-responsive gels Victor V. Yashin, Olga Kuksenok, Anna C. Balazs ∗ Chemical Engineering Department, University of Pittsburgh, Pittsburgh, PA, United States

a r t i c l e

i n f o

Article history: Received 2 September 2009 Received in revised form 19 October 2009 Accepted 22 October 2009 Available online 5 November 2009 Keywords: Self-oscillating gels Belousov–Zhabotinsky reaction Elastodynamics. Chemo-mechanical waves

a b s t r a c t Polymer gels undergoing the Belousov–Zhabotinsky (BZ) reaction are unique materials because the polymer network can undergo autonomous oscillations in the absence of external stimuli. We describe theoretical and computational approaches that allow us to simulate these BZ gels and highlight our recent findings that reveal the rich dynamical behavior exhibited by these distinctive materials. In particular, we found that the nature of the oscillatory behavior is highly dependent on the size of the sample. Thus, by tailoring the sample’s size, we can design gels that exhibit a uniform swelling and deswelling or display more complicated shape changes, which involve complex chemo-mechanical traveling waves. We then discuss a remarkable form of mechano-chemical transduction in these gels, where the application of an external force can drive a non-oscillatory system to exhibit sustained oscillations, and autonomous rotational motion. Due to this sensitivity to mechanical deformation, the BZ gels could potentially be used as a smart, responsive coating, which transmits a global signal to indicate that the system has been impacted. Finally, we focus on heterogeneous gels, which encompass regular arrangements of BZ patches within a nonreactive polymer matrix. By varying the placement of these BZ patches within the matrix, we could modify the functionality of the material or introduce multi-functional behavior within a single sample. Overall, the examples indicate that our methodology provides a powerful technique for probing the properties of such soft active materials. © 2009 Elsevier Ltd. All rights reserved.

Contents 1. 2.

3.

4.

Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.1. Kinetics of the BZ reaction in gels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2. Dynamics of the BZ gel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3. Model parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4. Formulation of the gel lattice spring model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Results and discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1. Chemo-mechanical transduction in 3D homogeneous BZ gels . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2. Mechano-chemical transduction in BZ gel films . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3. Modeling heterogeneous BZ gels: chemical communication in complex media . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Acknowledgement . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . References . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

∗ Corresponding author. E-mail address: [email protected] (A.C. Balazs). 0079-6700/$ – see front matter © 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.progpolymsci.2009.10.003

156 157 157 158 159 160 162 162 167 168 171 172 172

156

V.V. Yashin et al. / Progress in Polymer Science 35 (2010) 155–173

1. Introduction Over the last few years, there has been a surge of interest in the area of “soft active materials”, which transduce energy into mechanical work. 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 these networks can be driven to undergo significant swelling and deswelling. Typically, such rhythmic modulations are produced through controlled variations of the surrounding solvent [1–5] and the resulting oscillating gels can be used as micro-actuators [6] and pulsatile drug delivery devices [3,7]. There exists one system, however, where the polymer gel can expand and contract without external control or stimuli. This unique material is the so-called “BZ gel”, which comprises a chemo-responsive polymer network undergoing the Belousov–Zhabotinsky (BZ) reaction. Discovered in the 1950s, the BZ reaction is recognized as being a cornerstone of the field of nonlinear chemical dynamics [8,9]. The original reaction occurred in a fluid; later, chemically neutral, non-responsive polymer gels were used as a medium for the BZ reaction in order to suppress hydrodynamic effects [10]. It was only in the late 1990s that Yoshida et al. fabricated the first BZ reaction to take place in a responsive gel [11]. The BZ gels can undergo autonomous, self-oscillatory behavior due to the presence of a ruthenium catalyst, which is covalently bonded to the polymers [11–28]. The BZ reaction generates a periodic oxidation and reduction of the anchored metal ion. The hydrating effect of the oxidized metal catalyst induces an expansion of the gel, which then contracts when the catalyst is in the reduced state. In other words, the chemical energy from BZ reaction is tranduced into the mechanical oscillations of the gels. These self-oscillating BZ gels open up exciting opportunities for designing small-scale, biomimetic devices that operate in an essentially autonomous manner. In particular, such devices could perform sustained work until the reagents in the host solution are consumed; the system, however, can be readily “re-fueled” by simply replenishing these solutes. Millimeter-sized pieces of the BZ gels can actually undergo self-oscillations on the order of hours without replenishment of reagents [26,27]. Thus, the BZ gels could be harnessed to create novel micro-actuators, “soft robots”, and active components for microfluidic devices that exhibit self-sustained movement. In addition to their potential practical utility, the BZ gels provide an ideal medium for investigating non-linear dynamical phenomena that arise from a coupling of chemical and mechanical energy. To probe these phenomena and establish design rules for creating potentially useful systems, we recently undertook the first computational studies of BZ gels [29–37]. Modeling the behavior of such chemically responsive gels introduces a number of challenges. In particular, the chemical reaction alters the local polymer–solvent interactions; this leads to the redistribution of solvent inside the gel and gives rise to a deformation of the material. The latter deformation, in turn, affects the chemical reaction within the gel. Modeling this coupled

behavior is complicated by the fact that gel deformations are typically relatively large. To address the above issues, we recently formulated an approach that we refer to as the “gel lattice spring model” or gLSM, which allows us to model chemoresponsive gels undergoing finite deformations in two [30,31] or three [33] dimensions. The approach also allows us to uncover morphological transitions that arise from the interplay between chemical reactions and mechanical action [32,34,37]. By using the gLSM approach to simulate the BZ gels, we uncovered a rich variety of dynamic patterns and nonequilibrium behavior, and demonstrated that the method provides a powerful technique for probing the properties of such soft active materials [29–37]. In this review, we highlight some of our recent findings, focusing on examples that highlight the complexity of these chemo-responsive gels. In particular, we first describe how simply changing the size of a three-dimensional sample leads to dramatically different oscillatory behavior. Thus, by tailoring the size of the sample, we can design gels that exhibit a uniform swelling and deswelling or display more complicated shape changes, which involve complex traveling waves. We then discuss how the dynamical behavior of BZ gel films can be controlled by mechanically deforming the sample. Specifically, by applying an external force, we isolated cases where the polymer network can potentially act as artificial skin, which transmits a global signal to indicate that the system has been impacted. These studies on mechano-chemical transduction also revealed a scenario where the applied force caused the entire sample to undergo autonomous rotational and translational motion. Finally, we also describe our findings on the distinctive behavior of heterogeneous gels, which encompass regular arrangements of “BZ patches” in a non-reactive polymer matrix. Through these studies, we found that the functionality of the material can be altered by changing the arrangement of the BZ patches within the polymer matrix. In effect, the introduction of the BZ patches can be viewed as inscribing a “Lego set” within a polymer gel and using these inscribed components to build the desired device. In the ensuing discussion, we comment on the relevant experimental work where appropriate and refer the reader to a recent review by Yoshida [17] on current experimental efforts in creating novel BZ gel systems. While, to the best of our knowledge, there have been no prior computational studies on BZ chemo-responsive gels, there have been a number of theoretical models for probing the behavior of other types of oscillating gels [1,38–42], i.e., gels that oscillate in response to changes in the external environmental conditions. These studies provided significant insight into the factors that contribute to the regular pulsations. The prior calculations, however, were carried out in effectively one dimension since the systems were assumed to be spherically symmetric. While these 1D calculations can describe the volumetric changes in the system, they cannot account for changes in the shape of the sample. Thus, our 2D and 3D gLSM simulations provided the first approach for capturing the shape changes that can occur in oscillating gels and open up the possibility of uncovering new morphological transitions within the gels. Below, we

V.V. Yashin et al. / Progress in Polymer Science 35 (2010) 155–173

157

begin by providing a description of this gLSM model and then discuss how this general approach can be integrated with the governing equations for the BZ chemical reaction to simulate the behavior of BZ gels.

dY = −k1 H2 AY − k2 HXY + 12 fk5 BZ dt

(2)

dZ = 2k3 HAX − k5 BZ dt

(3)

2. Methodology

In a swollen polymer gel, the polymer network with the attached metal ion catalyst occupies a volume fraction . The rest of the BZ reagents are dissolved in solvent (of volume fraction 1 − ). The variable concentrations X, Y, and Z are defined with respect to the total volume of the system. In contrast, it is more convenient to define the concentrations of the reactants A, B, and H with respect to the volume of the solvent; in accordance with the original Oregonator model, we assume the concentrations of the latter reactants to be constant. Re-calculation of A, B, and H relative to the total volume is achieved through multiplication by the factor (1 − ). As a result, the following substitutions should be made in the right-hand-sides of the reaction rate equations Eqs. (1)–(3):

2.1. Kinetics of the BZ reaction in gels An example of a BZ gel is the cross-linked polymer network of poly(N-isopropylacrylamide) (NIPAAm) in which the catalyst ruthenium tris(2,2 -bipyridine) (Ru(bpy)3 ) is anchored onto the polymer chains [11–17]. This polymer network is swollen in an aqueous solution of NaBrO3 , HNO3 , and malonic acid (MA); MA acts as an organic substrate for the BZ reaction. While the kinetics of the BZ reaction involves two dozen variables for the concentrations of reactive species, as well as tens of chemical reactions, it is successfully described in terms of Field–Koros–Noyes (FKN) mechanism [43], which can be reduced to only three basic processes:

3

k1 H2 A → k1 (1 − ) H2 A

BrO3 − + 2Br− + 3H+ → 3HOBr

k2 H → k2 (1 − )H

BrO3 − + HBrO2 + 2Mred + 3H+ → 2HBrO2 + 2Mox + H2 O

k3 HA → k3 (1 − ) HA

2Mox + MA + BrMA → f Br− + 2Mred + other products

2

k5 B → k5 (1 − )B

The stoichiometric factor f represents the number of bromide ions produced as two oxidized metal ions (Mox ) are reduced (Mred ). In the reacting polymer gels [11–17], Mox = Ru(bpy)3 3+ and Mred = Ru(bpy)3 2+ . The cyclic change of the electric charge on the ruthenium ion causes the periodic swelling–deswelling of the gel in the course of the BZ reaction. The gel swelling takes place when the metal ion is in the oxidized state. The FKN mechanism is further approximated by the Oregonator model by the following sequence of reactions (the Tyson and Fife formulation [44]): k1 H2

A + Y −→X + P k2 H

X + Y−→2P

Thus, the polymer is assumed to act like a chemically neutral diluent, which affects the reaction rates only through changing the concentrations A, H, and B. In general, the stoichiometric factor f might also depend of the volume fraction . However, it is taken to be constant for simplicity. The concentration variable Y is known to change in time much faster than X and Z [45]. Consequently, the variable Y can be excluded from the reaction rate equations by employing the steady-state approximation dY/dt = 0 [44]. The remaining two reaction rate equations can be conveniently written in the dimensionless form after the dimensionless concentration variables u = X/X0 and v = Z/Z0 are introduced utilizing the parameterization given by Tyson [44,45]: X0 =

k3 HA , Z0 2k4

=

(k3 HA)2 . k4 k5 B

The dimension-

less time variable t = time/T0 is defined using the time scale −1 T0 = (k3 HA) . The resulting dimensionless reaction rate equations have the following form:

k3 H

A + X−→2X + 2Z k4

X + X−→A + P k5

B + Z−→ 12 f Y −

where A = [BrO3 ], B = [all oxidizable organic species], H = [H+ ], P = [HOBr], X = [HBrO2 ], Y = [Br− ], and Z = [Mox ]. The stoichiometric factor f is considered as a model parameter. The Oregonator treats the concentrations of the major reactants A, B, and P, as well as the concentration of the hydrogen ions H, as constants. The reaction rate equations for the species X, Y, and Z are as follows [44,45]: dX = k1 H2 AY − k2 HXY + k3 HAX − 2k4 X2 dt

(1)

du = F(u, v, ) dt

(4)

dv = εG(u, v, ) dt

(5)

where the reaction rate functions F and G from the Oregonator model derived by Tyson and Fife [44] are now modified to take into account the polymer volume fraction [29]: 2

F(u, v, ) = (1 − ) u − u2 − (1 − )f v 2

G(u, v, ) = (1 − ) u − (1 − )v

2

u − q(1 − )

2

u + q(1 − )

(6) (7)

158

V.V. Yashin et al. / Progress in Polymer Science 35 (2010) 155–173

The dimensionless parameters ε and q have the same definitions as in the original Oregonator model [44], 2k k k B namely: ε = k 5HA , q = k 1k 4 . 3

2 3

As can be seen from Eqs. (6) and (7), the Oregonator model specifies the kinetics of the BZ reaction through three parameters, namely q, ε, and f. The parameter q is a dimensionless combination of only the reaction rate constants, so it remains unchanged under variations of the chemical composition at a constant temperature. In contrast, the Oregonator parameters ε and f are sensitive to the chemical composition of the BZ substrate, and other experimental conditions. The parameter ε is a ratio of two reaction rates: one reaction rate is for reduction of the metal ion catalyst during the oxidation of the organic substrate and the other is for oxidation of the catalyst mediated by the BZ reactant u. Typically, oxidation of the BZ catalyst is faster than its reduction, so ε < 1. The stoichiometric parameter f specifies production of the bromide ion, which inhibits the catalytic process. Thus, the value of f effectively controls the concentration of oxidized catalyst within the system in the steady state and affects the amplitude of the oscillations in the oscillatory state [43]. It is common practice in modeling the BZ reaction to keep the value of q fixed, and consider the values of ε and f as the free model parameters, which control the reaction regime [43,45]. The modified Oregonator model, Eqs. (6) and (7), describes the kinetics of the BZ reaction in terms of the dimensionless concentrations of the oxidized catalyst, v, and the key reaction intermediate (the activator) u; the reaction rates explicitly depend on the volume fraction of polymer, . The spatial–temporal variations of  in the course of the swelling and deswelling of the gel affect the local chemical kinetics. Moreover, the BZ reaction itself induces the re-distribution of solvent with the swollen chemo-responsive hydrogel. In the next section, we will consider a theoretical model that couples the chemical reaction to the dynamics of the responsive polymer gel. 2.2. Dynamics of the BZ gel The BZ reaction in gel is accompanied by the diffusion of the dissolved activator u described by the diffusion flux j(u) , and by the transport of u and v with the solvent and polymer motions, which are characterized by the respective velocities v(s) and v(p) . Incorporation of these transport processes into the modified Oregonator, Eqs. (4) and (5), results in the following dimensionless reaction-diffusion equations: ∂u = −∇ · (uv(s) ) − ∇ · j(u) + F(u, v, ) ∂t

(8)

∂v = −∇ · (vv(p) ) + εG(u, v, ) ∂t

(9)

We note that Eq. (9) does not include a self-diffusion term because the catalyst is tethered to the chains. Eqs. (8) and (9) should be supplemented by the continuity equation for the volume fraction of polymer: ∂ = −∇ · (v(p) ) ∂t

(10)

The diffusion flux j(u) in Eq. (8) is calculated as that of a solute in a porous media with a porosity of (1 − ) [46]: −1

j(u) = −(1 − )∇ (u(1 − )

).

(11)

In this manner, the effect of the polymer on the diffusion of u is taken into account. The diffusion coefficient of the dissolved reagent, Du , is absent from the above equation because it is absorbed into the length scale L0 = (Du T0 )1/2 used to define the dimensionless coordinates in Eqs. (8)–(10). The velocity of the polymer network, v(p) , is determined by the local forces acting on the polymer chains. Due to a high polymer–solvent friction, the dynamics of the polymer network in swollen gels is often assumed to be purely relaxational (inertialess), so the forces that act on the deformed gel are balanced by the frictional drag due to the motion of the solvent [47,48]. The force balance equation can be written in the following dimensionless form:

∇ · ˆ = v0 T −1 Du ()(v(p) − v(s) )

[0, 1]

(12)

Here, ˆ is the dimensionless stress tensor measured in units of v−1 T , where v0 is the volume of a monomeric unit 0 within the polymer chains and T is temperature (measured in energy units and assumed to be constant throughout the system), and () is the friction coefficient. Both the polymer–solvent interactions and the elastic forces conˆ The friction coefficient is tribute to the stress tensor . taken to be of the form () ≈ 6s  −2 (), where the correlation length () scales like −3/4 [48]. This approximation for the friction coefficient is valid in the semidilute and intermediate regimes, i.e., for  < 0.5 [48], which is always the case in our calculations. To further simplify the model, we assume that only the polymer–solvent interdiffusion contributes to the gel dynamics and neglect the total velocity of the polymer–solvent system. Consequently, we specify that: v(p) + (1 − )v(s) = 0

(13)

If the stress tensor ˆ is known, then Eqs. (12) and (13) give the following equations for the polymer and solvent velocities to be used in the governing reaction-diffusion equations, Eqs. (8)–(10):

  −3/2

v(p) = 0 (1 − )

0

∇ · ˆ

(14)

−1 (p)

v(s) = −(1 − )

v

−1

Here, 0 = T (v0 (0 )Du ) is the dimensionless kinetic coefficient, and 0 is the volume fraction of polymer in the reference state, which is the undeformed state of the gel. The stress tensor  ˆ can be derived from the free energy density of the deformed gel, U, which consists of the elastic energy density associated with the deformations, Uel , and the polymer–solvent interaction energy density, UFH . ˆ and I3 = detB ˆ Namely, U = Uel (I1 , I3 ) + UFH (I3 ) where I1 = trB are the invariants of the left Cauchy-Green (Finger) strain ˆ [49]. The invariant I3 characterizes the volumetric tensor B changes in the deformed gel [49]. The local volume fractions of polymer in the deformed and undeformed states −1/2 are related in the following way:  = 0 I3 .

V.V. Yashin et al. / Progress in Polymer Science 35 (2010) 155–173

The elastic energy contribution Uel describes the rubber elasticity of the crosslinked polymer chains, and is proportional to the crosslink density, c0 , the number density of elastic strands in the undeformed polymer network. We use the Flory model [50] to specify Uel : Uel =

c0 v0 1/2 (I1 − 3 − lnI3 ) 2

(15)

It is worth noting that the polymer network in the BZ gels is sufficiently far away from the gel point (the typical value of Young’s modulus is approximately 2 × 104 Pa at 20 ◦ C (see ref. [51]) to justify the usage of the classical theory of rubber elasticity. The energy of the polymer–solvent interaction is taken to be of the following Flory-Huggins form [31]: UFH =

1/2 I3 [(1 − )

ln(1 − ) + FH ()(1 − )

∗

− v(1 − )],

(16)

where * > 0 describes the hydrating effect of the metal ion catalyst and captures the coupling between the gel dynamics and the BZ reaction, and FH () is the polymer–solvent 1/2

interaction parameter. In Eq. (16), the coefficient I3 appears in front of the conventional Flory-Huggins energy because the energy density is defined per unit volume of gel in the undeformed state. It is important to note that in Eq. (16), we neglect the electrolyte effects on the gel swelling behavior. This simplifying assumption is adopted because, according to the experimental study in ref. [26], the effect of variations in the ionic strength of the reacting BZ substrate on the gel swelling is much weaker than the hydrating effect of the oxidized metal ion catalyst. Finally, from Eqs. (15) and (16), one can derive the following constitutive equation for the chemo-responsive polymer gels [31]: ˆ = −P(, v)ˆI + c0 v0

 ˆ B, 0

(17)

where ˆI is the unit tensor, and the isotropic pressure P(, v) is defined as: P(, v)=−[+ln(1−)+ ()2 ]+c0 v0 (20 )

−1

+ ∗ v, (18)

where () = 0 + 1  is derived from the Flory-Huggins parameter, FH () [52]. We note that for non-responsive gel ( * = 0), the gel dynamics does not depend on the chemical reactions. In this case, our theoretical description of a swollen polymer gel is reduced to the model used by Hirotsu to analyze the experimental data in ref. [52]. 2.3. Model parameters To the best of our knowledge, a comprehensive characterization of a BZ gel that would be sufficient for indentifying all the above model parameters is currently absent in the literature. Recent review articles (refs. [16,17]) present results of experimental studies on BZresponsive gels, and provide a detailed description of the synthesis of the polymer network, the chemical composition of the reactive substrate, temperature, and other

159

Table 1 The values of the parameters of the 2+ Ru(bpy)3 -catalyzed Belousov–Zhabotinsky reaction at 20 ◦ C. BZ reaction in solutions, ref. [53] k1 = 2 M−3 s−1 k2 = 3 × 106 M−2 s−1 k3 = 42 M−2 s−1 k4 = 3 × 103 M−1 s−1 k5 = 6 M−1 s−1 BZ reaction in gels, ref. [12] A = 0.084 M B = 0.0625 M H = 0.3 M

experimental conditions. Data on the BZ reaction rate constants and on the polymer–solvent interaction, however, is lacking from the publications cited above. Nevertheless, a reasonable choice of model parameters to describe the BZ reaction in a hydrogel can be made based on data published on the Ru(bpy)3 2+ -catalyzed BZ reaction in solution, and on the swelling of NIPAAm gels in the absence of a chemical reaction [31]. Table 1 (see also ref. [31]) presents data on the reaction rate constants ki , 1 ≤ i ≤ 5, of the Oregonator model of the BZ reaction (see Eqs. (1)–(3)) published in ref. [53]. These values of ki were used to simulate the Ru(bpy)3 2+ -catalyzed BZ reaction in solution at the temperature of 20 ◦ C. This data can be utilized to model the BZ reaction in the polymer gel because the same catalyst is involved. The reaction kinetics also depends on the chemical composition of the BZ substrate through the concentrations A, B, and H (see Eqs. (1)–(3)). Table 1 presents the values of A, B, and H reported in the experimental study on traveling waves of swelling that were propagating along a specimen of the chemoresponsive gel [12]. Based on the data in Table 1, the values of the Oregonator parameters q and ε can be estimated as q = 9.52 × 10−5 and ε = 0.354. No experimental data is available to identify the value of the stoichiometric parameter f, so it is treated as an adjustable model parameter. Within the model of the chemo-responsive gel utilized here, the polymer backbone and the attached metal ion catalyst contribute to the polymer–solvent interactions through the interaction parameters () and * , respectively (see Eq. (18)). The parameter () can be specified based on the results of an experimental study on polyacrylamide gels without ionic or catalytic pendent groups [52]. It was demonstrated in ref. [52] that for a wide range of temperature, T, and polymer volume fraction, , the gel swelling could be described by taking the interaction parameter to be of the form (, T ) =

0 (T ) + 1 , where 0 is a certain function of the temperature T, and 1 is a constant. At the temperature of 20 ◦ C, the polymer–solvent interactions can be described by the function () = 0.338 + 0.518 [52]. The interaction parameter * , which mimics the hydrating effect of the oxidized metal ions, is an adjustable parameter of the model since no experimental data is currently available to guide our choice of * . The values of the volume fraction of the undeformed gel, 0 , and the dimensionless crosslink density in the gel, c0 v0 , in Eqs. ((17) and (18)) characterize the conditions of the

160

V.V. Yashin et al. / Progress in Polymer Science 35 (2010) 155–173

Fig. 1. Schematic of the 3D element. The entire sample consists of Lx × Ly × Lz nodes. In the underfomed state, the set of indexes i = 1 . . . Lx , j = 1 . . . Ly , and k = 1 . . . Lz defines the respective position of the nodes in x, y, and z directions. The node numbering within the element is marked by (1:8) and corresponds to the numbering within the entire sample (given here in parentheses) as follows: node 1 (i, j, k), node 2 (i, j, k + 1), node 3 (i + 1, j, k + 1), node 4 (i + 1, j, k), node 5 (i, j + 1, k), node 6 (i, j + 1, k + 1), node 7 (i + 1, j + 1, k + 1), node 8 (i + 1, j + 1, k). Coordinate system local to this element (,,) is marked by green arrows. Forces acting on node 1 (marked by the green circle) of element m = (i, j, k) are marked by the red and blue arrows. Namely, red arrows inside the element mark the spring-like elastic forces acting between node 1 and the next-nearest and next-next-nearest neighbors within the same element m (as defined in Eq. (22)). The blue arrows outside of the element mark contributions to nodal forces from the isotropic pressure within this element (as defined in Eq. (23)). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of the article.).

gel preparation. The experimental studies on the chemoresponsive gel reported in ref. [26] were conducted under the preparation conditions corresponding to 0 = 0.139 and c0 v0 = 1.3 × 10−3 . The transport processes in the gel are characterized by the diffusion coefficient, Du , of the dissolved reactant, u, in the solution, and by the dimensionless mobility coefficient of the gel, 0 . The value of Du = 2 × 10−5 cm2 s−1 can be used for the diffusion coefficient [43,45]. The numerical value of 0 is much greater that one at low values of the volume fraction of polymer [54]. We typically use 0 = 100 in the numerical simulations. Finally, we estimate the characteristic time and length scales of the problem. Under the chosen model parameters, −1 the time scale is T0 = (k3 HA) ∼1 s, and the length scale is 1/2 L0 = (Du T0 ) ∼40 ␮m. 2.4. Formulation of the gel lattice spring model To study the dynamic behavior of the BZ gels, we numerically integrate Eqs. (8)–(10) in two [30,31] or three [33] dimensions using our recently developed “gel lattice spring model” or gLSM. This method combines a finite element approach for the spatial discretization of the elastodynamic equations and a finite difference approximation for the reaction and diffusion terms. We used the gLSM approach to examine two-dimensional confined films and threedimensional bulk samples; here, we briefly discuss the more general 3D formulation [33]. We represent a 3D reactive, deformable gel by a set of general lineal hexahedral elements [55,56] (see Fig. 1). Initially, the sample is undeformed and consists of (Lx − 1) ×

(Ly − 1) × (Lz − 1) identical cubic elements; here Li is the number of nodes in the i-direction, i = x,y,z. The linear size of the elements in the undistorted state is given by . For a homogeneous BZ gel in the undeformed state, the polymer and crosslinks are uniformly distributed over the gel sample, i.e., each undeformed element is characterized by the same volume fraction 0 and crosslink density c0 . Upon deformation, the elements move together with the polymer network so that the amount of polymer and number of crosslinks within each hexahedral element remain equal to their initial values. Each element is labeled by the vector m = (i,j,k), and the element nodes are numbered by the index n = 1–8 and characterized by the coordinates rn (m) [33]. Within each element m, the concentrations of the dissolved reagent, u(m), the oxidized metal ioncatalyst, v(m), and the volume fraction of polymer, (m), are taken to be spatially uniform. The latter value is related to the volume of the element, V(m), as (m) = 3 0 /V (m). Since the polymer network dynamics is assumed to be purely relaxational, and Eq. (14) states that the network velocity v(p) is linearly proportional to the forces acting on ˆ the velocity of node n of the element m the network ∇ · , is proportional to the force acting on this node, Fn (m), i.e. drn (m) = Mn (m)Fn (m) dt

(19)

where Mn (m) is the mobility of the node. In refs. [31,33], we demonstrated that the nodal force can be readily determined by differentiating the finite element approximation of the energy of the deformed gel and the nodal mobility can be obtained by a proper discretization of the force balance equation, Eq. (12). Below, we briefly outline the major derivation steps, and refer the reader to the publications [31,33] for further details. In the finite element approximation, the total energy of the deformed system, Utot , is represented as a sum of the  elemental contributions, Utot = 3 m U(m), where U(m) is the elemental energy density defined per unit volume. Thus, the total energy can be expressed as a function of the coordinates of the nodes, and the force acting on a node can be calculated as Fn (m) = −∂Utot /∂rn (m). As a node can be common to several elements, the nodal force includes the contributions from all of the elements adjacent to the node. The elemental energy density U(m) depends on the deformations local to the element m. In general, the threedimensional deformation is specified by three functions: x = x(X,Y,Z), y = y(X,Y,Z), and z = z(X,Y,Z), where R = (X,Y,Z) and r = (x,y,z) are the respective coordinates of the same material point before and after deformation R → r(R). Inside each linear hexahedral element, these three functions are approximated by a tri-linear expansion. Within each element, we introduce a local coordinate system (, , ) as shown in Fig. 1. We then calculate the coordinates within the element m in this local coordinate system through the values of the nodal coordinates, rn (m), and a set of “shape functions”, as detailed in refs. [55,56]. We perform all the necessarily volume and surface integrations within each linear hexahedral element in this local coordinate system [33].

V.V. Yashin et al. / Progress in Polymer Science 35 (2010) 155–173

The energy density of the deformed gel U(I1 ,I3 ) (see discussion in Section II.B) can be expressed a as sum of 1/2 two terms, namely, U = u1 (I1 ) + u2 (I3 ). The first term describes a neo-Hookean elastic material, and is a linear function of the strain invariant I1 c0 v0 u1 (I1 ) = (20) (I1 − 3) 2 whereas the second contribution to U is a non-polynomial 1/2 function of J = I3 u2 (J) = UFH (J) −

c0 v0 lnJ. 2

(21)

The integration of the neo-Hookean energy density, Eq. (20), over the volume of element m can be carried out in the local coordinate system with no further approximations [33] to obtain the function U1 (m), which depends on the node coordinates rn (m). In contrast, Eq. (21) cannot be integrated explicitly, so we use the approximation U2 (m) ≈ u2 (−3 V (m)), where V(m) is the volume of element m as a function of rn (m) [33]. Finally, the total force acting  on a node is calculated by differentiating Utot = 3 m [U1 (m) + U2 (m)] with respect to the corresponding node coordinates. The total force acting on each node contains contributions from the elastic and osmotic properties of the system. It was shown [31] that the total force acting on node n of the element m consists of two contributions, namely, Fn (m) = F1,n (m) + F2,n (m). The first term, F1,n (m), originates from the neo-Hookean elasticity contribution to the energy, Eq. (20), and the second term, F2,n (m), originates from Eq. (21) and accounts for the isotropic pressure defined in Eq. (18). We have shown [33] that F1,n (m) can be expressed as a combination of linear spring-like forces, namely:

⎛

c0 v0  F1,n (m) = 12

+

Here,

⎝

 NN(m )





w(n , n)[rn (m ) − rn (m)]

and

NN(m )



(22)

represent the respective summa-

NNN(m )

tions over all the next-nearest neighbor nodal pairs and next-next-nearest neighbor nodal pairs belonging to all the neighboring elements m adjacent to node n of the element m. Above, w(n , n) = 2 in Eq. (22) if n and n belong to an internal face and w(n , n) = 1 in Eq. (22) if n and n belong to a boundary face [33]. Unlike the situation for purely two-dimensional deformations [31], there is no contribution from the interaction between nearest-neighbors in Eq. (22). The second contribution to the force acting on node n of the element m, F2,n (m) can be written as [33]: F2,n (m) =

1 P[(m ), v(m )][n1 (m )S1 (m ) 4

S(m)



Fn (m),

(24)

n

where n is the outward unit normal for an element’s surface, S(m). The total frictional drag force in the righthand-side of Eq. (12) can be calculated approximately by evaluating the integrand at the nodes of element m to obtain

dV0

  1/2 0

(1 − )

−1 (p)

v

3  n (m)vn (m). 8

(25)

n

Here, n (m) and vn (m) are the nodal friction coefficients and velocities, respectively, and we took into account the relationship between the deformed and non-deformed volume integration elements, namely, dV = 0 −1 dV0 . The frictional drag force depends on the local value of the volume fraction of polymer, . In the gLSM, the volume fraction of polymer is an attribute of an element as a whole. In order to evaluate the friction coefficient, n (m), at node n of element m, we average the volume fractions of polymer in the adjoining elements m and denote the result by

(m)n , i.e., for the internal node, (m)n = 18 m (m ). Finally, equating the right-hand-sides of Eqs. (24) and (25) yields the nodal mobility coefficient, Mn (m), in Eq. (19):

 (m) −1/2

Mn (m) = 8 0 −3 (1 − (m)n )

m

+ n2 (m )S2 (m ) + n3 (m )S3 (m )].

dSn · ˆ ≈

≈

NNN(m )





V0 (m)

⎞

[rn (m ) − rn (m)]⎠ .

m. The pressure within each element, P[(m ), v(m )], is calculated according to Eq. (18). In Eq. (23), the vector nl (m ) is the outward normal to the face l of element m , and Sl is the area of this face [33]. The vectors nl (m) are shown in Fig. 1 for element m that includes node n = 1 (here, faces 1, 2, and 3 correspond to  = −1,  = −1 and  = −1 in the local coordinate system, respectively). Both contributions to the force acting on node n = 1 of the element m from within this element are shown schematically in Fig. 1. The spring-like forces between node n = 1 and neighboring nodes are marked by red arrows, while the forces F2,n (m) are depicted by the blue arrows. We emphasize that the total force acting on node n of element m includes similar contributions from each of the neighboring elements containing this node. If the nodal forces are known, we can obtain the dynamics of the polymer network through Eq. (19). To obtain the mobility Mn (m), we integrate Eq. (12), which equates the network and the frictional drag forces, over the volume of element m, V(m) (here, we also take into account Eq. (13)). Integration of the left-hand-side of Eq. (12) gives the total force acting on an element m and is approximately equal to the sum of the nodal forces:

−1 0



161

(23)

In the above equation, the summation is performed over all the neighboring elements m’ that include node n of element

n

0

.

(26)

In order to obtain the equations for the concentration variables u(m) and v(m), we introduce the material time derivative dp /dt = ∂/∂t + v(p) · ∇ and re-write Eqs. (8) and

162

V.V. Yashin et al. / Progress in Polymer Science 35 (2010) 155–173

Fig. 2. Regular oscillations in chemo-responsive gel. (a–f) Snapshots of gel’s morphology during one period of oscillations at simulation times: (a) t = 1761, (b) t = 1767, (c) t = 1773, (d) t = 1779, (e) t = 1788, (f) t = 1794. The color represents the concentration of the oxidized catalyst, v, according to the color bar in (h). The minumum and maximum values for the color bar are vmin = 8 × 10−4 and vmax = 0.4166, respectively. (The same color bar (h) is used in all the following 3D images, whereas the values of vmin and vmax are given separatrely for each figure.) The size of the sample is 12 × 12 × 12 nodes and the stoichiometric factor f = 0.68. (g) Evolution in time of the average values < v > and  for the sample in (a–f); the points marked (a–f) correspond to the respective snapshots in (a–f).

(9) in the following form:



between the results obtained independently and those obtained with the gLSM.

u

dp u = −u∇ · v(p) + ∇ · v(p) dt 1−



3. Results and discussion



u + ∇ · (1 − )∇ + F(u, v, ) 1− dp v = −v∇ · v(p) + εG(u, v, ) dt

(27)

3.1. Chemo-mechanical transduction in 3D homogeneous BZ gels

(28)

In deriving Eq. (27), we took into account the equation for the diffusion flux of the reagent u, Eq. (11), and neglected the total velocity of the polymer–solvent system (see Eq. (13)). The local balance equations for u(m) and v(m) were specified by approximating the right-hand-sides of Eqs. (27) and (28) as described in ref. [33]. Namely, the chemical kinetics was described as being local to each element, the contributions from the polymer–solvent interdiffusion were calculated through the nodal velocities given by Eq. (19), and self-diffusion of the reagent u was approximated using the finite-difference technique. Ref. [33] provides a detailed description of the 3D numerical approach, and also contains a discussion of its validation. In particular, we considered a limiting case that can be solved via an independent method and demonstrated excellent agreement

As noted in Section 1, the unique feature of the BZ gels is their ability to autonomously convert chemical energy into mechanical action and it is in this sense that they effectively mimic biological systems. To develop a deeper understanding of this synthetic form of chemo-mechanical transduction, we simulated the threedimensional, dynamical behavior of the gels [33]. As described below, we found that the dynamics of this system depends not only on the reaction parameters, but also on the size of the sample. Herein, we present examples of this distinctive behavior and provide a phase map that pinpoints where the chemical oscillations produce regular or non-regular oscillations in the size and shape of the responsive gel. We also show that the dynamics of the BZ gels strongly depends on the boundary conditions at the surface of the sample [33].

V.V. Yashin et al. / Progress in Polymer Science 35 (2010) 155–173

163

Fig. 3. Effect of increasing the value of the stoichiometric factor in BZ reaction, f , on the dynamics of the sample. (a) Non-regular oscillations at f = 0.68; corresponding simulation times from top to bottom are: t = 1761, t = 1767, and t = 1770. The minumum and maximum values for the color bar (given in Fig. 2h) are vmin = 7 × 10−3 and vmax = 0.4342, respectively. (b) Regular oscillations at f = 0.9; corresponding simulation times from top to bottom are: t = 1761, t = 1788, t = 1797. The minumum and maximum values for the color bar are vmin = 3 × 10−4 and vmax = 0.2470, respectively. In both (a) and (b), the size of the sample is 24 × 24 × 24 nodes.

As the first example, we consider a BZ gel undergoing regular oscillations. The size of the sample is 12 × 12 × 12 nodes and we imposed the no-flux boundary conditions for u at the surface of the gel. (Unless specified otherwise, the ensuing examples also involve the no-flux boundary condition.) For the initial conditions, we chose the concentrations of the oxidized catalyst and reagent u to be randomly distributed around their stationary solutions, vst and ust , respectively. The gel can attain a steady state if the elastic stresses are balanced by the osmotic pressure and, simultaneously, the reaction exhibits a stationary regime for the same system parameters. More specifically, such stationary solutions (st , ust , vst ) can be found by numerically solving the following system of equations [33]: c0 v0

 1/3  st

0

st − 20

F(ust , vst , st ) = 0;



= osm (st , vst );

G(ust , vst , st ) = 0,

(29)

The left hand side of the first equation in Eq. (29) represents an elastic stress [31]. If the solution of the system of equa-

tions in Eq. (29) is known, the corresponding stationary 1/3 degree of swelling can be calculated as st = (0 /st ) . We used the latter value to set the initial size of each element, which was taken to be a cube with side st . Here,  = 1 in dimensionless units, which corresponds to L0 ∼ 40 ␮m in dimensional units, as described above. Fig. 2 shows the graphical output from our 3D gLSM simulations; the snapshots were taken during one period of oscillation at late times, ensuring that the simulations capture the regular, non-transient behavior. Within these images, the colors represent the concentration of oxidized catalyst, v, with the color bar being given in Fig. 2h. (We use the same color scheme in all the following 3D images, however, the values of vmin and vmax are different in each of the figures and are specified in the respective figure captions.) Fig. 2 also reveals the synchronization in time of the chemical and mechanical oscillations for this sample. As can be seen, the concentration of oxidized catalyst, v, attains its lowest value (see color bar) when the gel is in the least swollen state (Fig. 2a). As the concentration of oxidized catalyst increases in the course of BZ reaction,

164

V.V. Yashin et al. / Progress in Polymer Science 35 (2010) 155–173

the sample’s degree of swelling also increases, as shown in Fig. 2b–d. Correspondingly, when v begins to decrease, the degree of swelling also decreases, as illustrated in Fig. 2e–f. In Fig. 2g, we plot the evolution of v and , which are the respective average values of the concentration of the oxidized catalyst and the polymer volume fraction, for the sample shown in Fig. 2a–f. The average values are taken over all the elements within the sample at each moment of time. The dots marked (a)–(f) correspond to the images in Fig. 2a–f. The plot clearly shows that the lowest values of the average polymer volume fraction (red line) correspond to the highest values of the oxidized catalyst (black line), and visa versa. This behavior is similar to the in-phase synchronization of the chemical and mechanical oscillations observed experimentally by Yoshida et al. in cubic gel pieces that were smaller in size than the characteristic length scale of the chemical wave [15]. To gain further insight into the dependence of the dynamical behavior on the size of the sample, we next studied the dynamics of a sample of size 24 × 24 × 24 nodes; all other parameters were kept the same as in Fig. 2. The samples in Fig. 3 are larger then the characteristic length scale of the chemical wave, so that the effect of the diffusion of the reactive species within the sample is different here than in the example above. The snapshots in Fig. 3a show the sample’s evolution at late times and reveal that the oscillations are irregular; the actual realization of the dynamic pattern depends on the small random fluctuations in the initial conditions. What we observe is analogous to a fragment of a spiral wave that is typically observed in BZ reactive systems and the originating point for this wave depends on the small random fluctuations in the initial parameters of the system. If, however, we increase the value of the stoichiometric factor from f = 0.68 (in Figs. 2 and 3a) to f = 0.9, then the gel sample is observed to exhibit the regular periodic oscillations, as shown in Fig. 3b. The above examples illustrate that the behavior of the system strongly depends on both the size of the sample and the stoichiometric factor f. In the following simulations, we vary f from 0.4 to 0.95. The value of f is increased in increments of 0.05 except in regions close to the critical points, where we reduced these increments to 0.01 in order to more precisely define the critical values, f* and f** (see below). We conducted the simulations for four different sample sizes L × L × L (with L = 2, 6, 12 and 24 nodes). For each value of f and L, we ran three independent simulations with different random perturbations in the initial conditions. Each of these simulations was run for a sufficiently long time (until t = 2000) to ensure that the observed behavior is non-transient and robust for the given set of parameters. The findings from this series of simulations are summarized in the phase map in Fig. 4. The samples are observed to undergo regular oscillations for f > f* ; for these cases, we ˜ which is the average value of  around which plot , ˜ = (  the system oscillates in time. Specifically,  max +

min )/2, where max and min are the respective maximum and minimum values of  calculated at late times when the oscillations are regular [33]. On the other hand, the system is observed to be in the non-oscillatory regime when f ≤ fL∗ and for these systems, we plot the average

Fig. 4. Phase behavior for gel samples of various sizes. The sizes of the samples (L × L × L) nodes are provided in the legend; here, we choose L = 2, 6, 12 (solid lines) and L = 24 (diamonds). Black dashed line marked st represents a stationary solution (see Eq. (29)).

˜ that the sysvalue of the polymer volume fraction  =  tem reaches at late times. In the latter cases,  is equal to the value of the stationary solution, st , taken at the respective value of f and does not depend on the size of the sample. The black dashed line corresponding to st coincides with the lines connecting simulation data points for each of the values of L if f ≤ fL∗ (see Fig. 4). The plot in Fig. 4 clearly shows that increasing the size of the sample increases the critical value of fL∗ (for the range of sizes and the no-flux boundary conditions considered here). Thus, for example, a sample with L = 2 can exhibit oscillatory behavior at a lower value of f than a sample with L = 12. The latter fact has important implications if we consider a system that encompasses a number of gel samples of different sizes. In particular, if the stoichiometric factor ∗ , then all samples with L < 12 will undergo is set at f = f12 regular, non-decaying oscillations. On the other hand, all samples with L ≥ 12 will always reach the steady state at late times. Thus, by decreasing the size of the sample (while keeping the no-flux boundary conditions for u on the surface of the sample), one can induce transitions from the non-oscillatory to oscillatory regime [33]. Fig. 4 also reveals that for the samples with L ≤ 12, only one critical value, fL∗ , is observed in the simulations, such that for f > fL∗ , we observe regular oscillations of the sample. On the other hand, for larger-sized gels, the simulations yield two critical values. For example, when L = 24, we find that the system reaches a steady state at late times ∗ , and produces regular oscillations for f ≥ f ∗∗ ; for f ≤ f24 24 ∗ < f < f ∗∗ , the gel undergoes non-regular however, for f24 24 oscillations. An example of the latter behavior is shown in Fig. 3a. We mark the region of the non-regular oscillations for L = 24 with open circles in Fig. 4 because the ˜ (as defined above) only applies to regcalculation of  ular oscillations. Additional simulations showed that the region of non-regular oscillations also exists for samples of sizes L = 18 and L = 36. Our findings on the effect of the sample’s size on the dynamical behavior are similar to features observed

V.V. Yashin et al. / Progress in Polymer Science 35 (2010) 155–173

experimentally for spherical, non-responsive gel beads undergoing the BZ reaction [57]. In particular, the researchers observed [57] a switching between the “global oscillation” and “traveling wave” regimes of chemical oscillations as they increased the size of the spherical bead [57]. For the smaller bead sizes, the researchers observed uniform global oscillations, where the reactant concentrations at each moment of time were essentially uniform throughout the sample. This is similar to what we observe for the small sample sizes, where the contribution from diffusion is relatively small (i.e., for L ≤ 12). We note that for L = 12 the concentration of u is somewhat non-uniform as can be seen from Fig. 2; this non-uniformity, however, is rather small so that the effect of diffusion still remains small, i.e., in the above context, the oscillations can be regarded as “global”. With an increase in the size of the spherical bead, the same researchers observed a switching from the global oscillations to the traveling wave regime of chemical oscillations [57]. This behavior is similar to what we observe in our simulations for relatively large values of L (i.e., L = 24), as can be seen in Fig. 3. Finally, the authors also demonstrated [57] that the frequency of oscillations is higher for the cases of the traveling waves than for the cases of global oscillations. This again is qualitatively similar to our findings in ref. [33], where we showed that the period of oscillations is smaller for L = 24 (i.e., the case where we clearly observe traveling waves propagating throughout in the sample). We can make further qualitative comparisons between our findings and experimental results by investigating the effect of ε on the behavior of the BZ gels [13–15,25]. The value of ε is proportional to the concentration of malonic acid (MA) in the experiments. In the following simulations, we fix f = 0.8 and L = 12 and vary ε. The inset in the Fig. 5 illustrates the evolution of the average degree of

165

Fig. 5. Dependence of the period of oscillations on ε for the sample of size 12 × 12 × 12 nodes and with f = 0.8. Inset: evolution in time of the degree of swelling, (averaged over all the elements), for ε = 0.04 and ε = 0.6, respectively.

swelling in these samples for ε = 0.04 and ε = 0.6. The plot shows that the period of oscillations dramatically decreases with an increase in ε. This trend is also evident from the central plot, which shows the period of oscillations for a range of ε. The observed decrease in the period of oscillations with an increase in ε is in a qualitative agreement with the experimental results of Yoshida et al [13–15,25]. Finally, if we further increase the value of ε (i.e., ε ≥ 0.9 for the scenario presented above), we then observe a transition to the non-oscillatory regime. Again, this observation agrees qualitatively with experimental studies [13] where researchers observed a transition between the oscillatory and non-oscillatory regimes with an increase in the concentration of MA. Finally, we show that the dynamics of a gel sample depends strongly on the boundary conditions at the sur-

Fig. 6. Effect of the flux of u through the surface of the sample. (a and b) Evolution of the sample at the following simulation times: (a) t = 1767, (b) t = 1788. Images in (c) and (d) show only the bottom half of the sample in (a) and (b), respectively. Values for the color bar (Fig. 2h) are vmin = 6 × 10−5 and vmax = 0.3728.

166

V.V. Yashin et al. / Progress in Polymer Science 35 (2010) 155–173

Fig. 7. Compression-induced chemo-mechanical oscillations and rotations in 2D sample. (a) Schematic of the system. (b) Steady state for the sample with the degree of swelling in the vertical direction taken as ⊥ = 1.3. (c–f) Snapshots of the system after the compression to ⊥ = 1.18 at the following simulation times: (c) t = 8924, (d) t = 8927, (e) t = 8936, and (f) t = 8939 (here, we instantaneously decreased ⊥ from 1.3 to 1.18 at the simulation time t = 3 × 103 ). The size of this two-dimensional sample is 30 × 30 nodes, where the lattice unit in the unreformed state is taken to be  = 0.36; we also set 0 = 113 and f = 0.648. The color represents the value of v/ with (v/)min = 2 × 10−4 and (v/)max = 13.5324.

face of the sample. As we already noted, in all the three dimensional simulations presented above, we applied the no-flux boundary conditions for the dissolved reactant u at the surface of the samples [33]. Such no-flux boundary conditions correspond, for example, to the scenario where the entire sample is placed in a flexible casing, which does not permit a flux of reagent u in or out of the gel. In the next simulations, we consider another limiting case; namely, the concentration of u is kept at u = 0 outside of the sample and, therefore, there is a diffusive flux of u through the surface of sample into the outer solution. (To keep u = 0 in this outer solution, one should continuously remove the reactant u from the solution, as could be accomplished in continuously stirred tank reactors.) Fig. 6 shows snapshots of the late-time oscillatory behavior of BZ gels where we impose the u = 0 boundary conditions at the surface of the sample. The images in the top row show the distribution of v within the entire sample, and the images in the bottom row show the distribution of v in the bottom half of the sample only (i.e., the images in the bottom allow us to look “inside” the images in the top.) All the system parameters (including the initial conditions and the sample size) are identical to those in Fig. 3a; the only difference between the gel samples in Fig. 3a and Fig. 6 is the boundary conditions. By comparing Figs. 3a and 6, we observe that the sample with the no flux boundary condition undergoes nonregular oscillation and has a significantly higher degree of swelling than the sample with the diffusive flux through its boundary [33]. The diffusion of u through the surface of the

sample creates a gradient of u in the vicinity of the sample’s surface; this gradient leads to a decrease in the value of v in the course of the BZ reaction. This decrease in v induces the effective shrinking of the gel close to the surfaces in the sample. Hence, even though the initial size of the samples in Figs. 3a and 6 were chosen to be equal and the number of nodes was kept fixed at 24 × 24 × 24, the actual average size of the sample in Fig. 6 is smaller than in Fig. 3a, contributing to the fact that the oscillations of the sample become regular. The simulations discussed above illustrate that the type of oscillations observed in the system, as well as the period and the amplitude of these oscillations, depend strongly on the boundary conditions at the surface of the sample. This fact could potentially open up possibilities for controlling the oscillation of the gel by simply changing the concentration of u in the outside solution. In the above discussion, we primarily focused on the effect of the chemical oscillations on the mechanical pulsations of the gel. In the BZ gels, however, the transduction between chemical and mechanical energy is closely coupled, so that the gel’s mechanical oscillations also affect the chemical reaction. In particular, note that the evolution of u and v (Eqs. (8) and (9)) depend on , which varies as the sample undergoes mechanical swelling and deswelling. The BZ gels, however, can exhibit a more dramatic form of mechano-chemical transduction, where an applied mechanical force drives a stationary sample to exhibit pronounced oscillations. As a result, these BZ gels could potentially be used as sensors for mechanical impact,

V.V. Yashin et al. / Progress in Polymer Science 35 (2010) 155–173

167

Fig. 8. Non-decaying oscillations generated by local compression. Here, the sample size is 120 × 120 × 10 nodes, and f = 0.64. We apply an external force Fe = 1.6 × 10−2 to each element within an area of size of 12 × 12 elements. This force is applied along the “−z” direction as marked by the white arrow in (a). Corresponding simulation times are: (a) t = 1518, (b) t = 1740, (c) t = 1968, (d) t = 2202. The force was applied at t = 1.5 × 103 and removed at t = 2 × 103 , i.e., the snapshot in (d) is taken after the force is removed.

indicating both the magnitude and location of the applied force [37]. Below, we describe our findings on this distinctive interplay between mechanical action and chemical response. 3.2. Mechano-chemical transduction in BZ gel films In the next set of simulations, we focus on twodimensional gel samples and show that these films can be driven to exhibit dramatic and highly non-linear responses to relatively small compression. Fig. 7a shows a schematic of our system: a swollen, reactive gel is confined between the two parallel plates. Here, we assume the sample is uniform in the vertical z-direction, while it can freely swell and de-swell in the x- and y-directions. With these assumptions, the dynamics of the system can be regarded as purely two-dimensional; in the simulations, we use our 2D gLSM approach detailed in ref. [31]. The thickness of the sample in the z-direction is ⊥ H0 , where ⊥ is the degree of swelling in the z-direction, and H0 is the sample thickness in the undeformed state (see Fig. 7a). At a specified time, the slab is uniformly compressed and the degree of swelling in the z-direction is decreased by the relatively small amount, ␦ ⊥ ; we assume that this change in ⊥ occurs instantaneously [34]. Between such instantaneous compressions, the thickness of the sample remains constant. By considering confined BZ gels that are initially in the oscillatory state, we showed that a uniform compression dramatically affects the nature of the oscillations within the system [34]. We also showed that uniform compres-

sion can lead to even more pronounced effects if the sample is initially in a non-oscillatory steady state. Namely, for a selected parameter range, the compression of the sample can induce the chemo-mechanical oscillations, and such oscillations can, in turn, lead to the autonomous motion of the entire sample [34]. The images in Fig. 7 reveal the behavior of the latter case described above, where a uniform compression induces not only the oscillation of the gel, but also rotation of the entire sample [34]. In these simulations, we use the following parameters: f = 0.648, 0 = 113, and c0 = 1.3 × 10−3 ; the rest of the parameters are the same as presented above. Here, we also assume that u = 0 at the surface of the sample. We consider a square sample of size 30 × 30 lattice sites and the lattice unit in the undeformed state is taken to be  = 0.36 dimensionless length units. The initial system (with ⊥ = 1.3) exhibits the nonoscillatory steady state shown in Fig. 7b; the snapshots in Fig. 7c–f illustrate the subsequent compression-induced oscillations and rotation. These snapshots are taken during one half of the period of oscillations that occur after the sample was compressed to ⊥ = 1.18. As can be seen, a region with a high concentration of the oxidized catalyst moves counterclockwise within the sample. Such motion causes periodic shape changes and ultimately leads to the slow clockwise rotation of entire sample [34]. The spontaneous symmetry breaking in the compressed sample causes its rotation either clockwise or counter-clockwise; the direction that is actually chosen by the system depends on the small random fluctuations in the initial conditions. We note that the pressure required for compressing the

168

V.V. Yashin et al. / Progress in Polymer Science 35 (2010) 155–173

sample from ⊥ = 1.3 to ⊥ = 1.18 is estimated to be on the order of ∼4 × 103 Pa, which corresponds to a force of ∼6 × 10−3 N [34]. With respect to such autonomous motion, we isolated another example of BZ gel motility that is caused by the chemo-mechanical waves [31]. In particular, using the 2D gLSM simulations, we observed a BZ gel sample that migrated in the opposite direction from the propagating waves, which were generated by a pair of spirals [31]. In the above 2D systems, the swelling of the gel in the vertical direction is constrained to a constant value, so we could only model the effects of a spatially uniform deformation (e.g., uniform compression). Using our 3D simulations, however, we recently modeled the effect of applying a spatially localized force on a three-dimensional BZ film that is free to oscillate in the x, y and z directions [37]. Analogous to the case presented above, we found that a localized impact can drive a system that was initially in the nonoscillatory state into the oscillatory regime (see Fig. 8). The chemical waves are nucleated in the region of the local impact and propagate outwards. Within the 3D system, these variations in chemical concentration produce propagating “ripples” on the surface of the film [37]. Moreover, we isolated cases where the system continues to oscillate even after the release of the impact. Our results on these 3D systems provide the first predictions that local mechanical deformations can excite traveling chemical waves and wide-spread oscillations within BZ polymer gels [37]. The findings open up the possibility of harnessing BZ gels for a range of applications. Specifically, these materials could be used to create sensors that not only can transmit a signal in response to mechanical impact, but also transport reagents to address the after effects. Additionally, the BZ gels could be utilized to fabricate membranes that deliver reagents in response to an applied force; this feature provides a unique on–off mechanism for regulating the transport of compounds in the material. The coatings could potentially provide a “synthetic skin” for robotics. By probing the relationships between the features of mechanical stimuli and the activated chemical waves, we can facilitate the development of these devices. 3.3. Modeling heterogeneous BZ gels: chemical communication in complex media In the above studies, we assumed that the BZ gel encompassed a uniform distribution of catalyst throughout the system. Using our 2D gLSM approach, we introduced patches of BZ gels into a non-reactive polymer network. In other words, we assumed that the catalyst was localized in distinct domains. Our focus on such heterogeneous BZ gels is inspired by the structural heterogeneity and hierarchy found in responsive or adaptive biological materials. One example of this vital heterogeneity is provided by the structure of skeletal muscle, where the characteristic striations contribute to the functionality of this tissue. Another example is apparent in the morphology of bone, where the structural heterogeneity and hierarchy again contribute to the functioning of the material. By considering active gels that exhibit spatial heterogeneity, we can potentially estab-

lish guidelines for expanding the utility of these gels or integrating a number of functions (e.g., sensing, communication and actuation) into one sample. In our prior studies of heterogeneous BZ gels, we have in fact observed that the functionality of the material can be tailored by altering the location of the BZ patches [32]. As we show below, by simply moving the relative positions of the BZ patches within the non-reactive polymer, we can design either a miniature actuator, which displays concerted waves of expansion or contraction, or a pump, which drives fluids in a specific direction through the material. These findings point to a novel “modular” design approach that allows us to impart different functionality by varying the spatial arrangement of identical pieces of BZ gels within the matrix material. In carrying out these studies, our goal is to determine the optimal placement of the BZ patches to achieve the desired functionality. We note that our studies of heterogeneous gels are also motivated by an interest in understanding how chemical signals pass through complex media, which encompasses reactive and non-reactive domains [21,58–62]. In particular, how do the BZ pieces that are separated by non-reactive polymer communicate with each other? Through these studies, we can uncover the role that the non-reactive polymer network plays in the overall behavior of the material. As a step toward addressing these issues, we first considered a single BZ patch within a non-reactive matrix (see Fig. 9a). To clarify the terminology used here, we note that the BZ patch contains the anchored catalyst and the chemical oscillations only occur within this BZ patch. The bulk of the gel in the ensuing discussion (see Fig. 9) is non-reactive, i.e., no catalyst is anchored to the chains in this region. The BZ patch can be responsive or non-responsive. In the responsive region, * = / 0 and the gel responds to chemical reaction by undergoing mechanical oscillations. For the non-responsive case, * = 0 and the chemical reaction does not cause a swelling/deswelling of the gel, i.e., there is no chemo-mechanical coupling. The single BZ patch of responsive gel ( * = 0.105) in Fig. 9a is 5L0 in size and is surrounded by the non-reactive polymer network. To understand how the oscillations of the responsive domain affect the neighboring area, we focused on two elements, one at the center of the chemoresponsive BZ gel and another that is a distance of 5L0 to the right and within the non-reactive matrix. For each point, we calculated ␦ =  − eq as a function of time, where eq is the equilibrium volume fraction of polymer in the nonreactive gel. As anticipated, ␦ exhibits relatively large, periodic variations for the oscillating element (see Fig. 9b). The element in the non-reactive gel, however, also shows rhythmic changes in ␦, albeit smaller in magnitude than those for the reactive point. Thus, the pulsations in the reactive patch distort the neighboring matrix through the periodic uptake and release of solvent. The influence of the reactive domain on the surrounding matrix is also apparent in Fig. 9c, where we plot the concentration of u as a function of time for the two elements. Since the non-reactive matrix does not contain catalyst, there are no inherent chemical oscillations within this region. As can be seen in Fig. 9c, however, the test element within the nonreactive matrix does show rhythmic oscillations in u. The

V.V. Yashin et al. / Progress in Polymer Science 35 (2010) 155–173

Fig. 9. Chemo-mechanical oscillations at the center and outsize of a single BZ patch within a non-reactive matrix. (a) Schematic of the system; point 1 denotes the center of a square-shaped BZ patch of size 5L0 , and point 2 is located outside the patch at a distance of 5L0 from point 1. (b) Variations of  relative to eq at points 1 and 2. (c) Variations of u at points 1 and 2. The BZ patch is located within a square-shaped sample with a patch-toboundary distance of 10L0 . The gel dynamics was simulated at * = 0.105. The equilibrium volume fraction of polymer at * = 0 is eq ≈ 7.46 × 10−2 . The lattice spacing in the simulations is 1/2L0 . Here and in Figs. 10–13, we set f = 0.7, ⊥ = 1.1, and 0 = 100; the sample swells freely through the outer boundaries, and u = 0 outside the sample.

latter variations arise from the diffusion of u from the oscillating patch into the neighboring areas. (Note the relatively large changes in u within the BZ patch.) The plots in Fig. 9 reflect the fact that because the chains are cross-linked into a network, perturbations in one location are propagated through the system. The latter behavior will play an important role in facilitating communication between patches in systems that contain more than one reactive domain. We next consider the simplest multi-patch system, where two patches are separated by a non-reactive domain and the lateral separation between the patches is given by x (see Fig. 10). Both patches of responsive BZ gels exhibit periodic oscillations and, as can be understood from the above discussion, these oscillating patches interact, or effectively “communicate”, through the concentration

169

field u and the deformation of the gel around the reactive regions. Due to this interaction, the chemo-mechanical oscillations in the patches become synchronized [63]. We found that the two patches exhibit in-phase or out-of-phase synchronization. In a synchronized state, the phase difference between the two oscillators, , is locked to a specific constant value. For the in-phase synchronization (Fig. 10b),  = 0, while for the outof-phase mode (Fig. 10c) the phase difference is about  = . We now systematically varied the distance between the patches, x, while fixing * = 0.105 and found that the synchronization between the two oscillators depends critically on x. Fig. 11a shows that only the in-phase mode was found for x < 3, whereas both the in-phase and outof-phase modes were observed for x ≥ 3. Furthermore, the frequency of oscillation, ω, in a synchronized state also depends on x, and this dependence is remarkably different for the two modes. (Note that the frequency is normalized by ω0 , which is the frequency of oscillation within an isolated patch.) In particular, ω in the out-ofphase mode increases as the BZ patches are placed closer together until this mode ceases to exist at x < 3. In contrast, ω in the in-phase mode decreases with a decrease in x. Moreover, there exists a range of the inter-patch distances (3 ≤ x ≤ 6 in Fig. 11a), where the out-of-phase patches oscillate with a higher frequency than the patches exhibiting the in-phase mode. The above behavior is quite distinct from the non-responsive case (i.e., * = 0), where the polymer does not deform due to the chemical reaction and the BZ patches only interact through the chemical degrees of freedom. Fig. 11b shows that in this case, the oscillation death regime [64,65] was found to exist in the system at 1 ≤ x ≤ 3. The absence of the oscillation death at * = 0.105 allows us to conclude that the chemo-mechanical coupling has a stabilizing effect on the synchronized oscillations. The dependencies shown in Fig. 11a provide valuable guidelines for designing patterned BZ gels that exhibit the desired dynamic behavior. To illustrate the utility of these design rules, we consider five catalytic patches that are arranged in a linear array. First, we focus on situations where all the patches are placed equidistantly, with x = 2 and x = 4, as shown respectively in Fig. 12a and b. For the two-patch case at x = 2, only the in-phase synchronization is supported by the system (see Fig. 11a). The simulations reveal that the corresponding array of five patches also exhibits the in-phase synchronization (see Fig. 12d). (The pattern is slightly bent because the end patches interact only with one neighboring patch, and consequently, they exhibit a slightly higher frequency than the internal patches.) For x = 4, both the in- and out-of-phase modes of pair-wise synchronization exist, but the out-of-phase has a higher frequency than the in-phase mode (Fig. 11a). It is well known that in an ensemble of coupled chemical oscillators, the dynamical mode having the highest frequency always dominates [59,60,62,66]. It is evident from the dynamical pattern in Fig. 12e that at x = 4, the neighboring patches do indeed exhibit the out-of-phase synchronization.

170

V.V. Yashin et al. / Progress in Polymer Science 35 (2010) 155–173

Fig. 10. Synchronization modes in the patterned BZ gel. (a) Schematic of the system; the black squares denote the catalyst patches 1 and 2. (b) The in-phase synchronization. (c) The out-of-phase synchronization. The patch size is 5L0 and the patch-to-boundary distance is 10L0 . The spatiotemporal behavior of u, v, and  was obtained at x = 5L0 , * = 0.105, and shown along the crosssection a-b. The gel dynamics is shown during the period of time of 100T0 . The lattice spacing in the simulations is 1/2L0 .

Fig. 11. The frequency of oscillations, ω, of the in-phase and out-of-phase synchronization modes as a function of the inter-patch separation x. (a) Responsive gel. (b) Non-responsive gel. ω0 is the frequency of oscillation within an isolated patch. The shaded area in (b) indicates the domain of the oscillation death regime.

We can now utilize the above findings to design materials that exhibit a specific functionality. In particular, we now arranged the five catalyst patches of responsive BZ gel as shown in Fig. 12c; within this array, the distance between the two patches at the right end is x = 4, while all other reactive domains are spaced at x = 2. As anticipated from the findings in Fig. 11, the pair of patches at x = 4 oscillates out-of-phase, with a frequency that is greater than that of the in-phase mode at x = 2. The simulations show that this pair forms a pacemaker and thereby generates a chemical wave that propagates to the left (see Fig. 12f). This wave can be harnessed to create a micropump, which drives dissolved species from one end of the sample to another. The responsiveness of the gel to the BZ reaction is crucial for the existence of the dynamical behavior shown in Fig. 12d–f. A non-responsive polymer matrix ( * = 0) will not produce the rich dynamics discussed above. As seen in Fig. 12g–i, only the end patches exhibit the oscillations, whereas the oscillations within the three middle patches are quenched at * = 0. We also considered a horizontal BZ strip within a nonreactive gel matrix (see Fig. 13); in case 1, the strip is placed in the center of the sample and in case 2, this strip is placed at the edge (i.e., the left most, non-reactive gel is cut away). The dynamic patterns in the two samples are quite different. Specifically, in case 2, a wave is seen to propagate from the right to the left edge. The differences in the observed behavior can be attributed to the fact that the ends of strips experience different environments in the two scenarios. These exam-

V.V. Yashin et al. / Progress in Polymer Science 35 (2010) 155–173

171

Fig. 12. Dynamical patterns in the heterogeneous gel having five catalyst patches arranged in a linear array: the effect of the inter-patch spacing. (a–c) Schematic of the systems. (d–i) The spatiotemporal behavior of u along the central horizontal crosssection of the sample for (d–f) responsive and (g–i) non-responsive gels. The patch size is 5L0 and the patch-to-boundary distance is 10L0 . The lattice spacing in the simulations is L0 . The gel dynamics is shown during a period of time of 100T0 . Color bar: umin = 8.42 × 10−5 , umax = 0.357 at * = 0.105; umin = 8.22 × 10−5 , umax = 0.268 at * = 0.

Fig. 13. Dynamical patterns of the heterogeneous gel having a rectangular BZ patch: the effect of location of the patch within the sample. The patch is 33L0 long and 5L0 wide. The spatiotemporal behavior of u, v, and  is shown along the crosssection a-b. Case 1: Symmetric placement of the patch with the patch-to-boundary distance of 10L0 . Case 2: The left end of the BZ patch is placed at the sample edge, and the patch-to-boundary distance is 10L0 on the opposite side. The simulations were performed at * = 0.105. The lattice spacing in the simulations is L0 . The gel dynamics is shown during a period of time of 100T0 .

ples clearly reveal that the placement of a BZ patch within the sample plays an important role and can be used as a design tool. In particular, case 2 can be harnessed to create a pump that transports fluid and reagents to the edge of the gel. The above results also highlight the proposed modular design concept, where starting with identical BZ patches, one can introduce different functionality by simply varying the arrangement of the patches.

4. Conclusions Polymer gels undergoing the BZ reaction provide an ideal medium for investigating the rich dynamical behavior that can occur from the coupling of non-linear chemical reactions and mechanical deformations. In particular, we can probe chemo-mechanical transduction and the feedback that the mechanical response produces on

172

V.V. Yashin et al. / Progress in Polymer Science 35 (2010) 155–173

the chemical reaction–or mechano-chemical transduction. Furthermore, these materials allow us to examine a particularly fascinating form of mechano-chemical transduction, where an externally applied pressure causes the gel to produce traveling chemical waves. The behavior of the BZ gels yields significant insight into pathways that lead to the efficient inter-conversion of chemical and mechanical energy. Such insight can facilitate the creation of novel biomimetic materials that autonomously convert chemical energy into mechanical action and thereby perform a self-sustained function. For example, this form of chemo-mechanical transduction can be harnessed to create self-propelled soft robots and micromachines [18,28,67]. Our findings on mechano-chemical transduction can provide guidelines for designing sensors that transmit a signal in response to mechanical impact and transport reagents to address the after effects. The findings will provide the design rules for fabricating coatings for damage detection, and for membranes for selective transport. In particular, polymer gels provide optimal materials for membranes. By isolating the conditions where the system undergoes mechano-chemical transduction, one can design membranes that deliver reagents in response to an applied force. This feature provides a unique on–off mechanism for regulating the flow of compounds in a material. In the above, we also described the first studies to determine how the compartmentalization of BZ gels in a non-reactive polymer matrix affects the dynamical behavior of the system. The results point to a novel “modular” design approach, which can impart the desired functionality to the material. In particular, identical pieces of BZ gel form the crucial components and it is the arrangement of these pieces that determines and provides the desired performance. For example, the out-of-phase pattern shown in Fig. 12e can be utilized to encrypt a chemical message or code in the material. In contrast, the unidirectional wave motion shown in Fig. 12f can be used to pump reagents or fluid through the system, while the concerted behavior seen in Fig. 12d might be valuable for fabricating a microactuator. This approach can be viewed as inscribing a Lego set within a polymer gel and using these inscribed components to build the desired device. The computational studies described herein provide the necessary guidelines for assembling the components into the appropriate pattern to yield the specified function. Acknowledgement ACB gratefully acknowledges financial support from ARO. References [1] Pojman JA, Tran-Cong-Miyata Q. Nonlinear dynamics in polymeric systems. ACS Symp Ser, vol. 869. Washington, DC: ACS; 2004. [2] Shibayama M, Tanaka T. Volume phase-transition and related phenomena of polymer gels. Adv Polym Sci 1993;109:1–62. [3] Dhanarajan AP, Misra GP, Siegel RA. Autonomous chemomechanical oscillations in a hydrogel/enzyme system driven by glucose. J Phys Chem A 2002;106:8835–8. [4] Szanto TG, Rabai G. pH oscillations in the BrO3 − –SO3 2− /HSO3 − reaction in a CSTR. J Phys Chem A 2005;109:5398–402.

[5] Labrot V, De Kepper P, Boissonade J, Szalai I, Gauffre F. Wave patterns driven by chemomechanical instabilities in responsive gels. J Phys Chem B 2005;109:21476–80. [6] Howse JR, Topham P, Crook CJ, Gleeson AJ, Bras W, Jones RAL, et al. Reciprocating power generation in a chemically driven synthetic muscle. Nano Lett 2006;6:73–7. [7] Mujumdar SK, Bhalla AS, Siegel RA. Novel hydrogels for rhythmic pulsatile drug delivery. Macromol Symp 2007;254:338–44. [8] Belousov BP. Collection of short papers on radiation medicine. Moskow: Medgiz; 1959. [9] Zaikin AN, Zhabotinsky AM. Concentration wave propagation in two-dimensional liquid-phase self-oscillating system. Nature 1970:535–7. [10] Yamaguchi T, Kuhnert L, Nagy-Ungvarai Z, Mueller SC, Hess B. Gel systems for the Belousov–Zhabotinskii reaction. J Phys Chem 1991;95:5831–7. [11] Yoshida R, Takahashi T, Yamaguchi T, Ichijo H. Self-oscillating gel. J Am Chem Soc 1996;118:5134–5. [12] Yoshida R, Kokufuta E, Yamaguchi T. Beating polymer gels coupled with a nonlinear chemical reaction. Chaos 1999;9:260–6. [13] Yoshida R, Onodera S, Yamaguchi T, Kokufuta E. Aspects of the Belousov–Zhabotinsky reaction in polymer gels. J Phys Chem A 1999;103:8573–8. [14] Sakai T, Yoshida R. Self-oscillating nanogel particles. Langmuir 2004;20:1036–8. [15] Yoshida R, Tanaka M, Onodera S, Yamaguchi T, Kokufuta E. In-phase synchronization of chemical and mechanical oscillations in selfoscillating gels. J Phys Chem A 2000;104:7549–55. [16] Yoshida R. Design of functional polymer gels and their application to biomimetic materials. Curr Org Chem 2005;9:1617–41. [17] Yoshida R. Self-oscillating polymer and gels as novel biomimetic materials. Bull Chem Soc Jpn 2008;81:676–88. [18] Murase Y, Maeda S, Hashimoto S, Yoshida R. Design of a mass transport surface utilizing peristaltic motion of a self-oscillating gel. Langmuir 2009;25:483–9. [19] Shinohara S, Seki T, Sakai T, Yoshida R, Takeoka Y. Photoregulated wormlike motion of a gel. Angew Chem Int Ed 2008;47:9039– 43. [20] Shen J, Pullela S, Marquez M, Cheng ZD. Ternary phase diagram for the Belousov–Zhabotinsky reaction-induced mechanical oscillation of intelligent PNIPAM colloids. J Phys Chem A 2007;111:12081–5. [21] Tateyama S, Shibuta Y, Yoshida R. Direction control of chemical wave propagation in self-oscillating gel array. J Phys Chem B 2008;112:1777–82. [22] Maeda S, Hara Y, Yoshida R, Hashimoto S. Peristaltic motion of polymer gels. Angew Chem Int Ed 2008;47:6690–3. [23] Maeda S, Hara Y, Yoshida R, Hashimoto S. Control of the dynamic motion of a gel actuator driven by the Belousov–Zhabotinsky reaction. Macromol Rapid Commun 2008;29:401–5. [24] Suzuki D, Yoshida R. Temporal control of self-oscillation for microgels by cross-linking network structure. Macromolecules 2008;41:5830–8. [25] Suzuki D, Yoshida R. Effect of initial substrate concentration of the Belousov–Zhabotinsky reaction on self-oscillation for microgel system. J Phys Chem B 2008;112:12618–24. [26] Sasaki S, Koga S, Yoshida R, Yamaguchi T. Mechanical oscillation coupled with the Belousov–Zhabotinsky reaction in gel. Langmuir 2003;19:5595–600. [27] Miyakawa K, Sakamoto F, Yoshida R, Kokufuta E, Yamaguchi T. Chemical waves in self-oscillating gels. Phys Rev E 2000;62:793–8. [28] Maeda S, Hara Y, Sakai T, Yoshida R, Hashimoto S. Self-walking gel. Adv Mater 2007;19:3480–4. [29] Yashin VV, Balazs AC. Modeling polymer gels exhibiting selfoscillations due to the Belousov–Zhabotinsky reaction. Macromolecules 2006;39:2024–6. [30] Yashin VV, Balazs AC. Pattern formation and shape changes in selfoscillating polymer gels. Science 2006;314:798–801. [31] Yashin VV, Balazs AC. Theoretical and computational modeling of self-oscillating polymer gels. J Chem Phys 2007;126:124707. [32] Yashin VV, Balazs AC. Chemomechanical synchronization in heterogeneous self-oscillating gels. Phys Rev E 2008;77, 046210/1–7. [33] Kuksenok O, Yashin VV, Balazs AC. Three-dimensional model for chemoresponsive polymer gels undergoing the Belousov–Zhabotinsky reaction. Phys Rev E 2008;78, 041406/1–16. [34] Kuksenok O, Yashin VV, Balazs AC. Mechanically induced chemical oscillations and motion in responsive gels. Soft Matter 2007;3:1138–44. [35] Yashin VV, Van Vliet KJ, Balazs AC. Controlling chemical oscillations in heterogeneous Belousov–Zhabotinsky gels via mechanical strain. Phys Rev E 2009;79, 046214/1–12.

V.V. Yashin et al. / Progress in Polymer Science 35 (2010) 155–173 [36] Dayal P, Kuksenok O, Balazs AC. Using light to guide the self-sustained motion of active gels. Langmuir 2009;25:4298–301. [37] Kuksenok O, Yashin VV, Balazs AC. Global signaling of localized impact in chemo-responsive gels. Soft Matter 2009;5:1835– 9. [38] Yui N, Mrsny R, Park K, editors. Reflexive polymers and hydrogels: understanding and designing fast responsive polymeric systems. Boca Raton, FL: CRC Press; 2004. [39] Boissonade J. Self-oscillations in chemoresponsive gels: a theoretical approach. Chaos 2005;15, 023703/1–9. [40] Boissonade J. Simple chemomechanical process for self-generation of rhythms and forms. Phys Rev Lett 2003;90, 188302/1–4. [41] Boissonade J, De Kepper P, Gauffre F, Szalai I. Spatial bistability: A source of complex dynamics. From spatiotemporal reactiondiffusion patterns to chemomechanical structures Chaos 2006;16, 037110/1–16. [42] Boissonade J. Oscillatory dynamics induced in polyelectrolyte gels by a non-oscillatory reaction: a model. Euro Phys J E 2009;28:337–46, doi:10.1140/epje/i2008-10425-1. [43] Scott SK. Oscillations, waves and chaos in chemical kinetics. New York: Oxford University Press; 1994. [44] Tyson JJ, Fife PC. Target patterns in a realistic model of the Belousov–Zhabotinskii reaction. J Chem Phys 1980;73:2224– 37. [45] Tyson JJ. A quantitative account of oscillations, bistability, and traveling waves in the Belousov–Zhabotinskii reaction. In: Field RJ, Burger M, editors. Oscillations and traveling waves in chemical systems. New York: Wiley; 1985. p. 93–144. [46] Chadam J, Peirce A, Ortoleva P. Stability of reactive flows in porous media: coupled porosity and viscosity changes SIAM. J Appl Math 1991;51:684–92. [47] Onuki A. Theory of phase-transition in polymer gels. Adv Polym Sci 1993;109:63–121. [48] Barriere B, Leibler L. Kinetics of solvent absorption and permeation through a highly swellable elastomeric network. J Polym Sci Part B Polym Phys 2003;41:166–82. [49] Atkin RJ, Fox N. An introduction to the theory of elasticity. New York: Longman; 1980. [50] Hill TL. An introduction to statistical thermodynamics. Reading, MA: Addison-Weley; 1960. [51] Maeda S, Hara Y, Yoshida R, Hashimoto S. A chemo-mechanical rotational actuator driven by BZ reaction. IEEE 2006:1160–5.

173

[52] Hirotsu S. Softening of bulk modulus and negative Poisson’s ratio near the volume phase transition of polymer gels. J Chem Phys 1991;94:3949–57. [53] Amemiya T, Ohmori T, Yamaguchi T. An Oregonator-class 2+ model for photoinduced behavior in the Ru(bpy)3 -catalyzed Belousov–Zhabotinsky reaction. J Phys Chem A 2000;104: 336–44. [54] Achilleos EC, Christodoulou KN, Kevrekidis IG. A transport model for swelling of polyelectrolyte gels in simple and complex geometries. Comp Theor Polym Sci 2001;11:63–80. [55] Smith IM, Griffiths DV. Programming the finite element method. Chichester, UK: Wiley; 2004. [56] Zienkiewicz OC, Taylor RL. The finite element method, vol. 1. Oxford, UK: Butterworth-Heinemann; 2000. [57] Aihara R, Yoshikawa K. Size-dependent switching of the spatiotemporal structure between a traveling wave and global rhythm. J Phys Chem A 2001;105:8445–8. [58] Steinbock O, Kettunen P, Showalter K. Anisotropy and spiral organizing centers in patterned excitable media. Science 1995;269:1857–60. [59] Mikhailov AS, Showalter K. Control of waves, patterns and turbulence in chemical systems. Phys Rep-Rev Sec Phys Lett 2006;425:79–194. [60] Fukuda H, Tamari N, Morimura H, Kai S. Entrainment in a chemical oscillator chain with a pacemaker. J Phys Chem A 2005;109:11250–4. [61] Taylor AF, Tinsley MR, Wang F, Huang ZY, Showalter K. Dynamical quorum sensing and synchronization in large populations of chemical oscillators. Science 2009;323:614–7. [62] Kheowan OU, Mihaliuk E, Blasius B, Sendina-Nadal I, Showalter K. Wave mediated synchronization of nonuniform oscillatory media. Phys Rev Lett 2007;98, 074101/1–4. [63] Pikovsky A, Rosenblum M, Kurths J. Synchronization. A universal concept in nonlinear sciences. Cambridge, UK: Cambridge University Press; 2001. [64] Crowley MF, Field RJ. Electrically coupled Belousov–Zhabotinskii oscillators. 1. Experiments and simulations. J Phys Chem 1986;90:1907–15. [65] Crowley MF, Epstein IR. Experimental and theoretical studies of a coupled chemical oscillator: phase death, multistability and in-phase and out-of-phase entrainment. J Phys Chem 1989;93:2496–502. [66] Kuramoto Y. Chemical oscillations, wave, and turbulence. Berlin: Springler-Verlag; 1984. [67] Maeda S, Hara Y, Yoshida R, Hashimoto S. Self-oscillating gel actuator for chemical robotics. Adv Robot 2008;22:1329–42.