Chemomechanical origin of directed locomotion driven by internal chemical signals

See allHide authors and affiliations

Science Advances  01 May 2020:
Vol. 6, no. 18, eaaz9125
DOI: 10.1126/sciadv.aaz9125


Asymmetry in the interaction between an individual and its environment is generally considered essential for the directional properties of active matter, but can directional locomotions and their transitions be generated only from intrinsic chemical dynamics and its modulation? Here, we examine this question by simulating the locomotion of a bioinspired active gel in a homogeneous environment. We find that autonomous directional locomotion emerges in the absence of asymmetric interaction with the environment and that a transition between modes of gel locomotion can be induced by adjusting the spatially uniform intensity of illumination or certain kinetic and mechanical system parameters. The internal wave dynamics and its structural modulation act as the impetus for signal-driven active locomotion in a manner similar to the way in which an animal’s locomotion is generated via driving by nerve pulses. Our results may have implications for the development of soft robots and biomimetic materials.


Diverse, often complex, locomotion is a key feature of active matter far from equilibrium (14). Elements of matter consume energy to interact with other elements and/or with the environment. Global locomotion and other collective behaviors then organically emerge on a macroscopic level. For living organisms, locomotion typically involves (5) fundamental direct and retrograde waves, i.e., motion propelled by muscular waves along or opposed, respectively, to the direction of motion, and more sophisticated modes such as walking, swimming, reciprocating migration, etc. A general understanding of the dynamic origin of active matter locomotion has not yet been achieved (1). A solution to this challenging task would be a major contribution toward the design of soft robots and active materials.

Asymmetry in the interaction between an individual and its environment is generally considered essential for the directional properties of active matter locomotion (14, 6). This principle also applies to the locomotion of living organisms; e.g., gastropods (such as snails) interact with substrates through viscoelastic mucus (with hysteretic asymmetry properties), while annelids (e.g., earthworms) harness their setae to produce an asymmetric frictional force during peristaltic interactions with the ground (5). In artificial systems, three methods are frequently used to introduce asymmetries into the system. In the first, asymmetries are incorporated into individual units (e.g., Janus particles) and generate (7) self-diffusiophoresis, self-electrophoresis, self-thermophoresis, self-acoustophoresis, or self-osmophoresis to drive motion. In the second, uniform particles or samples are subject to an inhomogeneous external field, such as a gradient of electric field, temperature, solute concentration, acoustic pressure, etc., that generates stresses and flows near the particle surface to drive its movement (8). Photoresponsive materials in the presence of patterned light (9) are another example of this approach. Moreover, the coupling dynamics achieved by using a field gradient can produce autonomous locomotion. Examples include the self-shadow effect of a photoactive polymer film to drive autonomous peristalsis of the film under directional illumination (10), the positive and negative phototactic locomotion of a Belousov-Zhabotinsky (BZ) gel under differential illumination (11), and the periodic migration of a BZ gel attributed to autonomous transitions of chemical dynamics under linear gradient illumination (12). The third method is a combination of the above two approaches, i.e., coupling between a heterogeneous field and a nonuniform medium. An example is the use of a time-varying magnetic field to drive and control a small untethered sample containing an embedded magnetic structure to perform the complex, multimodal locomotion involved in swimming, climbing, rolling, walking, jumping, or crawling (13). Note that although locomotion of the magnetic sample exhibits complex functionality, it is not autonomous.

Nevertheless, there are many examples in nature indicating that the origin of active matter locomotion is not necessarily an asymmetric interaction between an individual and its environment or an evolving asymmetrically embedded structure. For example, an animal’s locomotion can be not only complicated by limbs or wings, but it is also related to the system dynamics of internal interactions that generate stress asymmetries and play the basic role of a driving force (5). To cite two examples, amoeboid locomotion is caused by signal communication between individual cells (14), and the autonomous movement of a liquid crystal droplet is due to the collective dynamics of active nematics inside the drop (15).

These points raise the following fundamental questions: In addition to exploiting external asymmetries in natural and man-made systems, can locomotion and transitions between locomotive modes be generated solely by the intrinsic dynamics of the system (16)? Can this general problem concerning active matter locomotion be studied in a relatively simple chemical system, regardless of the much greater complexity of living organisms and fluid systems?

Here, we try to answer these questions by studying a bioinspired active gel (17) locomoting in a homogeneous environment. Our system is a widely studied self-oscillating polymer gel that hosts the tris(bipyridine)ruthenium [Ru(bipy)]-catalyzed Belousov-Zhabotinsky reaction (BZR). It is a reaction-diffusion-stress coupling system, which can transform chemical energy into mechanical work autonomously, producing not only periodic swelling-deswelling but also continuous locomotion of the gel (18). The redox state of the catalyst undergoes periodic changes during the homogeneous BZR oscillations, causing the gel to undergo autonomous swelling and deswelling, reminiscent of the rhythmic muscle contractions in animals (17). In a large enough gel sample, the propagation of chemical waves (spatiotemporal patterns of the concentration of oxidized catalyst) drives gel peristalsis, resulting in locomotion of the whole sample. This phenomenon is analogous to nerve pulses driving the locomotion of animals through deformable muscles (19). The gel lattice-spring model (gLSM) proposed by Balazs and coworkers (2022) captures the chemomechanical dynamics of the BZ gel, i.e., the large-scale shape changes (23), locomotion, and collective behaviors of the gel sample (24).

The first effort toward understanding BZ gel locomotion reported (25) that the direction of motion results from a counter-effect of polymer-solvent interdiffusion during the chemical waves, which propels the gel against the direction of wave motion (retrograde wave locomotion). Subsequently, under a step-changed light intensity condition, another fundamental mode of locomotion, i.e., direct wave locomotion, and a retrograde-direct wave transition, were found, which was attributed to a “chemical anchoring effect” modulating the push-pull competition of the chemical waves (19). That study revealed that the internal dynamics (push-pull effects) of the system play a dominant role in the gel locomotion. A more complex form of locomotion, reciprocating migration of the gel, was subsequently found (12), in which an autonomous transition between retrograde and direct wave locomotion is rooted in the interaction between the system’s dynamic instability and a gentle spatial gradient of illumination. Like neuronal signals, which can be described, for example, by the Fitzhugh-Nagumo (2628) or Hodgkin Huxley model (29), BZ waves, described by the Oregonator model, interacting with an active gel are nonlinear signal waves, which can serve as an analog of active (e.g., animal) locomotion.

It is, however, difficult to extract the essential origin of the gel locomotion from these earlier studies (11, 12, 19) due to the presence of differential illumination (i.e., the gel interacts with a heterogeneous environment), which introduces asymmetrical interaction into the system. In this work, we address the questions raised above by simulating a BZ gel that can move in a periodic space under a homogeneous environment; that is, no asymmetrical effect is introduced into the system from external conditions. We also examine the effects of modifying the parameters that govern the internal dynamics of the gel.


Locomotion modes and their transition in one spatial dimension

To investigate the gel locomotion in a homogeneous environment, we use pulse waves to propel a stimulus-responsive gel under uniform light intensity on a pseudo–one-dimensional (1D) ring. The 1D gel system can be seen as a collective of homogeneous synchronizing lattice oscillators. Once initialized, e.g., by a random fluctuation, the phase differences between these diffusion mechanics–coupled oscillators evolve to generate pulse waves that repeatedly propagate through the periodic boundary even in the absence of illumination, and no asymmetry is introduced into the gel sample via the boundary conditions or the environment (see fig. S1). The coupled oscillators lead to propagating pulse waves, resulting in directional locomotion of the entire gel. To speed up the simulations, we initiate pulse waves to propagate in the gel with the most likely wave number, i.e., one-pulse waves, according to fig. S2. A modified gLSM is used to study the gel locomotion, in which the kinetics of the photosensitive BZR, i.e., photoinduction and photoinhibition effects (30), the relationship between extreme values of the v variable and the light intensity (19), and the illumination-induced bifurcation of local oscillations, are described by a two-variable photosensitive Oregonator model (11, 31). The chemical kinetics and the structure of the chemical waves can be modulated by varying the intensity of the homogeneous illumination (I), which allows us to control the gel locomotion propelled by the chemical waves. The details of the model, parameters, calculation procedures, and simulation settings are given in Materials and Methods.

We next consider the chemical dynamics and locomotion of a photosensitive BZ gel modulated by increasing levels of homogeneous illumination. Key transitions occur at the following light intensities Ia = 0.0625, Ib = 0.0809, Ic = 0.0870, and Id = 0.117: (i) simple-to-complex pulse waves with interval-asymmetry structures (IASs), (ii) retrograde structured pulse waves to a stationary balance between the driving effects of primary and backfired waves, (iii) gel with no net motion to direct wave motion, and (iv) breaking up of the IASs. These phenomena are discussed in more detail below.

At low light intensity (0.0 ≤ I < Ia), e.g., I = 0.02, as shown in Fig. 1A, simple pulse waves repeatedly propagate through the gel, propelling the gel to move opposite to the waves (inset view of the net displacement shown in Fig. 1B), generating retrograde wave gel locomotion. One can identify two wave regions in these simple waves: the wave front and the wave back are indicated in Fig. 1A. These regions respectively push and pull the local gel to move along the +x and −x directions alternately (“1-push” and “2-pull” in the inset of Fig. 1B), and the net displacement over 1 cycle is in the −x direction, as shown in the inset of Fig. 1B. We emphasize that gel locomotion is driven by simple pulse waves even without illumination, i.e., I = 0, indicating that internal interaction between BZ waves and active matter is the origin of the stress asymmetry between push and pull.

Fig. 1 Spatiotemporal plots of v and the motion of the gel center under homogeneous illumination.

The two columns display spatiotemporal plots of v and position of the gel center versus time, respectively. The light intensity is (A and B) I = 0.020, (C and D) I = 0.070, (E and F) I = 0.087, and (G and H) I = 0.110. Color denotes the concentration of v. Labels 1 to 4 in (G) denote the approximate positions of the four wave regions in one interval structure, respectively, the primary wave’s front and back and the backfired wave’s front and back. Movie S1 shows the gel locomotion at I = 0.0, 0.02, 0.070, 0.087, and 0.116.

When the intensity exceeds Ia, IAS begin to emerge in the simple pulse waves. As shown in Fig. 1, C, E, and G, the simple pulse wave is transformed into a complex pulse wave. An IAS is a structural unit of these complex pulse waves produced by an illumination-induced spatiotemporal distortion. A complete IAS is shown in the inset to Fig. 1G. Each complex pulse wave is formed by the arrangement of six IASs (NIAS = 6 under illumination intensities Ia < I < Id), and each IAS is made up of two local waves, that is, a primary wave and a backfired wave. Primary waves are local waves that propagate along the overall direction of the pulse waves (red arrow in Fig. 1G, inset), whereas backfired waves are local waves that propagate in the opposite direction (blue arrow in Fig. 1G, inset). Each IAS has four wave regions, the primary wave’s front and back and the backfired wave’s front and back, which are roughly located at numbers 1 to 4, respectively, in Fig. 1G. With increasing illumination (Ia < I < Id), each IAS undergoes gradual distortion, as shown in Fig. 1 (C, E, and G). Note that the overall direction of pulse wave propagation does not change during this process. Under the driving force of these gradually distorting waves, the dynamics of the locus of the gel’s center, Rc, (adjusted for the moving baseline) changes gradually from quasiperiodic (Fig. 1D) to aperiodic oscillation (Fig. 1, F and H) with increasing illumination.

Meanwhile, the net displacement of the Rc locus shows the gel locomotion transition from a retrograde wave (I = 0.070, Fig. 1D) through a motionless state (I = Ic, Fig. 1F) and, ultimately, to a direct wave (I = 0.110, Fig. 1H). This transition is mediated by the IAS of the waves, and the zero net-displacement state of the gel is the critical state of the transition (light intensity Ic is the transition point). In this work, we have not found direct waves driven by simple pulse waves, which indicates that the presence of IAS may be necessary for the transition between retrograde and direct wave locomotion. When the intensity of illumination exceeds Id, the IAS breaks up into reverse-propagating pulse waves (see fig. S3), and the gel returns to retrograde wave locomotion comparable to that seen at low I, but in the opposite direction. Note that the retrograde-direct transition can also be generated by adjusting other system parameters such as ε (Oregonator kinetic parameter), f (Oregonator stoichiometric parameter), χ* (chemistry-gel dynamic coupling parameter), K (wave number) and Lgel (size of the gel), as shown in figs. S4 to S7.

Kinematic analysis of gel locomotion

To clarify the origin of the homogeneous light-induced transition between locomotion modes, we quantify key characteristics of the pulse wave and the gel locomotion by introducing several parameters: the local curvature [κ(x)] of the chemical waves (12), the mean velocity of the gel movement (Vg), and the mean net displacement of the gel center for each wave region (Li). The first step needed to ensure that these characterizations are useful is the division of the pulse wave into regions. We accomplish this by using the total differential of v in the spatiotemporal plots, which captures the integrated dynamics of the pulse waves. The differential expression readsdv=vxdx+vtdt(1)

Using Eq. 1, we propose a criterion to divide the pulse wave regions, as shown in Table 1.

Table 1 A criterion for identifying the subregions of chemical waves.

View this table:

We choose this criterion based on the following reasoning. On the one hand, at the wave fronts of both simple (Fig. 1A) and complex (Fig. 1, C, E and G) pulse waves, the local value of v increases with time (vt > 0), while at all wave backs, it decreases with time (vt < 0). On the other hand, when the overall motion of the pulse waves lies along the +x direction (with I < 0.117), the spatial gradients of v (i.e., vx) are negative and positive, respectively, at the primary wave front and wave back, but for backfired waves, the signs of vx are reversed due to the wave traveling in the −x direction, as shown in Fig. 1G. Therefore, the wave’s regions can be numerically identified by the criterion for vt and vx in Table 1. Using this criterion, the changes of the wave regions with increasing I may be distinguished, as shown in Fig. 2A. Simple pulse waves (e.g., the plot with I = 0.02) contain only the primary wave’s front (1) and back (2). For complex pulse waves (plots with I = 0.07 to 0.116), the backfired wave front (3) and wave back (4) are adjacent to the primary wave front and wave back.

Fig. 2 Pulse wave–driven gel kinematics under different intensities of homogeneous illumination.

Spatiotemporal plots of (A) v and (B) Vp. Red and black dashed lines, respectively, denote vt = 0 and vx = 0, which coincide for simple pulse waves, as shown in the first column of spatiotemporal plots of (A) and (B). (C) Spatial distribution of κ(x) versus I. The plot consists of 116 data points at an interval of ΔI = 0.001. (D) Vg versus I. (E) Lprimary and Lbackfire versus I. (F) Li versus I, i = 1 to 4. Light intensities Ia to Id denote critical points as follows: appearance of backfired waves, Lprimary = Lbackfire, Vg = 0, and IAS breakup, respectively. I, II, and III denote regions of simple, complex, and multiple waves, respectively.

For gel locomotion driven by pulse waves, the mean velocity of gel movement (Vg) represents the push-pull effect of the waves with increasing light intensity. Integrating the instantaneous velocity of the local gel, i.e., Vp, as shown in Fig. 2B, over each wave region yields the displacement in that region, LiLi=1Sj=1Sk=1G(Vp)iΔt(2)where S denotes the grid number of the gel, G = Tt, T is the period of the pulse wave, and (Vp)i denotes velocity of grids located in the wave region i. The sum of Li over all IAS in one wave is the net displacement (LT) driven by one pulse signalLT=i=1NLi(3)

Then, Vg is given byVgT=LT(4)where N denotes the number of characteristic regions of a pulse wave (for a simple pulse wave, N = 2, and for a complex wave, N = 4).

Using the above characteristic parameters, we analyze the transition from retrograde to direct wave locomotion of the gel. Note that the growing instability of the simple pulse waves with increasing illumination is the key process modulating the driving force acting on the gel. Therefore, we first examine how the structure of the chemical waves varies with the light intensity. For chemical waves, a spatial plot of the local curvature [κ(x)] versus I describes the distortion of the waves with increasing illumination, as shown in Fig. 2C. The defining equation for κ(x) can be found in Eq. 11 of Materials and Methods. In Fig. 2C, when I < Ia (region I), κ(x) is always zero, which indicates that only simple pulse waves propagate in the gel. As I increases from Ia to Id (region II), a spatially periodic organization of κ(x) begins to emerge and gradually strengthen, leading to six periodic structures (IAS) and distortion of the pulse waves. Meanwhile, comparing Fig. 2A with Fig. 2C, we see that when I exceeds Ia, wave regions 3 and 4 start to appear. The shape of the waves shortly after onset is seen at I = 0.07, near the transition point (Ia), in the spatiotemporal plots shown in Fig. 2 (A and B). The plots show a backfired wave front (wave region 3), which stems from the primary wave front (wave region 1). Its back (wave region 4) originates from the primary wave’s back (wave region 2). The complex dynamics in wave region 2 (primary wave back) play a key role in the transition of the locomotion mode and will naturally arise in our subsequent analysis.

The partial displacements Li (i = 1, 2, 3, 4) serve as the blocks from which the gel locomotion is built, since they characterize the propulsive contribution from each wave region. The key lies in finding the spatiotemporal region that dominates the mode transition of the gel locomotion. Fig. 2D, the Vg-I curve, shows the net velocity of gel locomotion as a function of the illumination intensity in simple pulse wave–driven retrograde wave locomotion [region I, Vg < 0 and κ(x) = 0] as well as in complex pulse wave–propelled retrograde wave locomotion [region II, Vg < 0 and κ(x) ≠ 0] and complex pulse wave–driven direct wave locomotion [region II, Vg > 0 and κ(x) ≠ 0]. From Eqs. 3 and 4, and the above discussion, we haveVgT=Lprimary+Lbackfire(5)where Lprimary = L1 + L2 and Lbackfire = L3 + L4 are the net displacements driven by the primary waves and the backfired waves, respectively. For a simple pulse wave, there is no Lbackfire. To examine the competition between the primary and backfired waves that determines the transition from negative to positive Vg, we plot Lprimary and Lbackfire versus I in Fig. 2E. Here, the light intensity Ib denotes the critical point where Lprimary = Lbackfire. Near the retrograde-direct transition point Ic, Lprimary increases roughly linearly with I, while Lbackfire (which is negative) decreases at a rate that diminishes with I. A clearer view is shown in the inset of Fig. 2E, where |Lprimary|/|Lbackfire| increases linearly across Ic and is always larger than 1, indicating that Lprimary dominates the transition from negative to positive Vg at Ic. The decomposition of Lprimary and Lbackfire into L1 to L4 with increasing illumination around Ic is shown in Fig. 2F. One can see that L1 is nearly constant, whereas |L2| decreases substantially, resulting in a negative-to-positive transition of LT, indicating that L2 dominates the change of Lprimary. On the other hand, for backfired waves, the two curves L3-I and L4-I are roughly symmetric about the zero line (horizontal dashed line in Fig. 2F), so that they almost cancel each other. In summary, wave region 2, which is the primary wave back, is the key region for the retrograde-direct transition of the gel locomotion.

In this stimulus-responsive gel with chemomechanical activity, the locomotion of the gel stems from the chemically induced asymmetry of the internal stress over the sample. We now analyze the internal dynamics of Vp by examining the effect of homogeneous illumination on the kinematic terms, i.e., the mean values of the gel velocity [(Vp)meani], the stress gradient [(σx)meani], and the apparent friction coefficient [(ξ′)meani] in each wave region, as shown in Fig. 3. The relationship between the local velocity of the gel and the stress is given by (21)Vp=σxζ(ϕ)(6)where σx=(πosmσel1D), ζ(ϕ)=Λ01(1ϕ)(ϕ/ϕ0)3/2, Vp is the grid velocity of the gel, σx is the stress gradient, and ζ′(ϕ) is the apparent friction coefficient. The mean values of the kinematic parameters [(Vp)meani] (σx)meani, and (ζ)meani] represent the average over each wave region (Fig. 3, A, C, and D). Ai, the area that each wave region occupies in the plots (Fig. 3, A and B), is also shown in Fig. 3B.

Fig. 3 Effect of light intensity on the kinematic parameters in each region of an IAS.

Mean values of kinematic parameters in each wave region: (A) (Vp)meani versus I. (B) Ai versus I. (C) (σx)meani versus I. (D) (ξ′)meani versus I.

Each displacement Li can be decomposed into (Vp)meani × Ai/S, as shown in Fig. 3 (A and B). (Vp)meani varies linearly with illumination near Ic in all four wave regions. However, (Vp)mean2 has a markedly smaller slope. The areas of the wave regions also undergo a roughly linear change with I around Ic; i.e., A3 and A4 increase, while A1 and A2 decrease. Our analysis of (Vp)meani and Ai suggests that the kinematics of L2 are dominated by the velocity of the gel rather than the relative areas of the wave regions.

We next analyze the mechanics of the local gel velocity (Vp) under the driving of the chemical waves using Eq. 6. The mean value of the apparent friction viscosity (ζ′)meani is investigated in Fig. 3D and fig. S8G. Under illumination Ic (see fig. S8G), the distribution of (ζ′)meani over each wave region is relatively uniform. Changes in the sign of Vp arise from σx rather than from ζ′, which is always positive. We can estimate (Vp)meani in each region by averaging Eq. 6(Vp)imean(σx)imean(ζ(ϕ))imean(i=1,2,3,4)(7)

As shown in Fig. 3C, each curve of (σx)meani versus I shows a trend that parallels the corresponding (Vp)meani-I curve in Fig. 3A. Each (ξ′)meani changes little with I, especially for the primary wave. As the inset to Fig. 3D shows that the ratio of (ξ′)mean1 to (ξ′)mean2 is almost independent of I, implying that the notable decrease of |L2| near Ic (see Fig. 2F) is dominated by the mechanical term (σx)mean2. Note that positive and negative σx values favor direct and retrograde waves, respectively; thus, the changes in the pull effect [(σx)mean2] with illumination facilitate the transition from retrograde to direct waves.

The above analysis suggests that the insensitivity of (σx)2mean to changes in I is the key to the locomotion transition. As shown in fig. S8 (E and F) (at I = Ic), Vp and σx display areas of both push and pull behavior within wave region 2. It seems clear that the Vp-induced decrease of |L2| (see Fig. 2F) and the σx-induced decrease of |Vp| in wave region 2 (Fig. 3A) arise from the same mechanism. Both are generated from the internal competition between opposing domains in wave region 2. In Fig. 4A, at the critical light intensity (Ic), we can distinguish three sub-areas (2a, 2b, and 2c) in wave region 2 in which σx is either negative (2a) or positive (2b and 2c).

Fig. 4 Mixed push and pull effects within primary wave back and dynamic bifurcation diagram of maximum of v at the gel center.

(A) Sub-areas of σx plot at I = Ic. Red and black dashed lines, respectively, denote vt = 0 and vx = 0. (B) Spatial profiles of v, ϕx, vx, and vxx within wave region 2 at I = Ic along the vertical red line in (A). (C) Dynamic bifurcation diagram of maximum of v at the gel center (vmax-center) versus I. The insets are the trajectories of the u-v oscillations at different light intensities.

As the illumination level increases, the net force that generates the retrograde-direct wave transition arises from two effects. First, there are four subregions, that is, wave regions 1, 2a, 2b, and 2c, in the primary wave of complex pulses. Both push (wave region 1) and pull effects (wave region 2a) in the primary wave are strengthened with illumination (see fig. S8F). Second, with increasing illumination, two additional subregions (2b and 2c) contribute growing push effects to offset the pull effect from subregion 2a, leading to a nearly constant value of (σx)mean2 with increasing illumination. Briefly, the mixing of the push and pull effects within wave region 2 leads to weakening of the total pull effect and, ultimately, to transition from retrograde to direct wave gel locomotion.

Dynamic origin of the transition from retrograde to direct wave locomotion

Why do multiple areas of push and pull effects coexist within wave region 2? We seek here to answer this question by examining the dynamical details of the stress gradient. In the gLSM, σx is described by the following equationσx=B(ϕ,v)ϕxχ*ϕvx(8)whereB(ϕ,v)=11/(1ϕ)+2χ0ϕ+3χ1ϕ2χ*v+C0(λϕ0/ϕ21/(2ϕ0))

The right side of Eq. 8 contains two composite terms, B,vx and ϕvx (χ* is set to 1.0), where ϕ and v are dynamical variables of the system. The quantities, v, ϕ, vx, ϕx, and B are shown for pulse waves at several illumination intensities in fig. S8 (A to D and H, respectively). Sub-areas of wave region 2 are marked as 2a, 2b, and 2c in each of these plots. From fig. S8H, we see that B is always negative. The volume fraction, ϕ, is positive by definition. As shown in Fig. 4B and fig. S8C, vx is always positive in wave region 2. So, according to Eq. 8, for σx to be positive, ϕx must be negative, i.e., we require a pushing force. Further, Fig. 4B shows that negative ϕx and small positive vx occur in areas 2b and 2c of the pulse waves to give a push effect, offsetting the pull effect in area 2a of the pulse waves.

The dynamic term ϕx in region 2 of the pulse wave is the key factor for the mode transition of the gel locomotion. As shown in Fig. 4B, ϕx is antiphase to vx. The evolutions of ϕ, ϕx v, and vx obey reaction-diffusion-mechanics coupling dynamics according to Eq. 9. As shown in Fig. 4C, the dynamics of local oscillations at the gel center undergo a bifurcation as the level of illumination is increased, changing from simple oscillations (e.g., the inset showing the u-v limit cycle with a stable periodic orbit at I = 0.020) to complex oscillations (e.g., the u-v trajectories at I = 0.070, 0.087, and 0.110), with the bifurcation occurring at Ib. Because the pulse waves are generated from the diffusion of local oscillations through the gel, such a complex local oscillation appears to be necessary to support the backfired waves (see fig. S8A). In addition, as shown in the insets of Fig. 4C, when I > Ib, with increasing illumination the orbit expands in the u-v plane, which leads to a greater variation of v over time and space. The changes in ϕ are strongly correlated with those in v, which accounts for the stimulus-responsive nature of the gel. Thus, over time, an illumination-enhanced oscillatory amplitude of v leads to an increasing amplitude of ϕ, ϕx, and vx oscillations (see fig. S9, C and D). In essence, the parameter-induced instability of simple local oscillations constitutes the dynamic origin of the backfired waves, which leads to the transition from retrograde to direct wave locomotion of the gel in the absence of external asymmetries.


In summary, we used a modified gLSM model to study the locomotion of a mechanosensitive BZ oscillating gel in a periodic space. Two locomotion modes, that is, retrograde and direct wave locomotion, as well as the transition between them, can be found in such a symmetrical system by adjusting kinetic or mechanical parameters. The numerical results indicate that the retrograde-direct transition of the gel locomotion occurs only if the gel is driven by complex pulse waves with IAS. Detailed analysis reveals that the complex pulse waves host four wave regions, including two push and two pull areas. Changes in the structure of these regions as the parameters are varied modulate the competition between the push and pull effects to generate the transition between modes. We have found that the dominant factor lies in the behavior of the primary wave back, in which homogeneous light-enhanced structural changes of the IAS inhibit the pull effect from this wave region to produce the retrograde-direct transition. Analysis of the chemical dynamics uncovers the origin of the locomotion, which is that the simple pulse waves lose stability to a more complex waveform, as modulated by the homogeneous illumination intensity and/or other parameters.

Further, in the absence of external asymmetric environmental interaction, or even in an isolated system without illumination (movie S1), autonomous active locomotion can be driven by interaction between simple internal nonlinear pulses and active responsive matter within the system. This interaction can be modulated by system parameters to change the locomotion direction and velocity, i.e., interaction with an external asymmetric environment is not a prerequisite for active locomotion. Here, by adjusting the system parameters, complex waves (local complex oscillations or chaos) result in reversal of the direction of locomotion. It is unlikely that such a transition can be realized without the loss of stability of simple pulse waves to more complex waves.

Coupling between nonlinear pulse signals, responsive active matter and a structured boundary [e.g., asymmetric system construction (7, 10) or appendages (32) like setae, mucus, limbs, or wings] should result in more diverse and sophisticated autonomous locomotion and collective behavior without the need for external stimulation.

Last, adaptability and “intelligent” motion of active matter in response to environmental stimuli have been reported, such as dynamic asymmetric light (9) and magnetic field–induced (13) versatile locomotion. Addition of internal nonlinear driving signals and active population communication would increase such systems’ computational ability to respond to environmental stimuli, producing functionally complex and potentially evolutionary motion that would enhance the level of intelligence in active matter. One can envision that designed DNA signals with different chemical reaction networks (33, 34) might be coupled with mechanical forces to construct a biocompatible active system for delivering and releasing drugs by using both internal and external control.


Our BZ gel model is constructed by coupling the Yashin-Balazs model (gLSM) (2022) with the two-variable form of the Amemiya model of the photosensitive BZR (11, 31). The governing partial differential equations (PDEs) for a 1D gel are{dϕ/dt=ϕVpdu/dt=uVp+(uVp/(1ϕ))+((1ϕ)(u/(1ϕ)))+F(u,v,ϕ,I)dv/dt=vVp+G(u,v,ϕ,I)(9)Vp=σ/ζ(ϕ)σ=(ϕ+ln(1ϕ)+χ(ϕ)ϕ2χ*ϕv+c0v0[(ϕ/ϕ0)1λ4ϕ(2ϕ0)1])ζ(ϕ)=(Λ0(1ϕ))1(ϕ/ϕ0)3/2where (11, 2022, 31) u and v are the nondimensionalized concentrations of HBrO2 and Ru(bpy)33+ in the Oregonator model, respectively; ϕ is the volume fraction of gel; ∇σ denotes the stress gradient, c0 is the cross-link density, and v0 denotes the volume of a monomeric unit within the polymer chains, with c0v0 = 1.3 × 10−3. Vp is the velocity vector of the polymer. The apparent friction coefficient, ζ′(ϕ), is a function of ϕ. Λ0 is a dimensionless kinetic coefficient set to 100. χ* is a coupling constant that couples the chemical reaction to the gel dynamics. ϕ0 denotes the undeformed polymer volume fraction, which was set to 0.139. The functions F and G describe the reaction kinetics of the BZR (modified Oregonator model). We replace the simple Oregonator in the Yashin-Balazs model with a two-variable photosensitive Oregonator model that incorporates both photoinduction and photoinhibition in the Ru(bipy)-catalyzed BZR. The functions F(u,v,ϕ,I) and G(u,v,ϕ,I) that describe the gel-coupled BZR are{F(u,v,ϕ,I)=(1ϕ)2uu2(fv+IP1)(1ϕ)[uq(1ϕ)2[u+q(1ϕ)2]1+IP2G(u,v,ϕ,I)=ε[(1ϕ)2u(1ϕ)v+I(0.5 P1+P2)](10)

Here, ε, f, and q are the Oregonator parameters. I denotes the dimensionless homogeneous light intensity. P1 and P2 are rate coefficients that characterize the two key photochemical reaction steps. Our simulation of the gel uses 200 grid points (denoted by S), and the gel length (Lgel) is set at 126.0 U. All times, space, variables, and parameters are dimensionless (7, 12, 21). The choice of model parameters ensures an oscillatory gel medium.

The PDEs were numerically integrated using the gLSM computational approach (2022), using an unequal distance grid approximation for the 1D Laplacian operators (22) (first and second derivatives). The integration time step was Δt = 1.0 × 10−4, and the space step was varied according to the local λ//, where λ// characterizes the degree of swelling along the unrestricted direction of the gel along the long(x) axis of the tube. The other two (restricted) dimensions are characterized by λ, which is fixed. Initially, the BZ gel was assumed to be in the steady state. The initial volume fraction was ϕ1D ini = ϕ0/(λ//λ2). The default values for our model simulation were f = 1.0, ε = 0.3, χ* = 1.0, P1 = 0.124, P2 = 0.77, λ// = 1.26, and λ = 1.06.

The waves can arise spontaneously due to random fluctuations. Initially, the variables u, v, and ϕ were assigned their stationary values (uss, vss, and ϕss) (21, 22), where uss and vss are the solutions (22) of F(uss, vss, ϕss) = 0, G(uss, vss, ϕss) = 0 at χ* = 0.0. We then added random noise (50% of the amplitude of uss) to all grid points, which induces synchronous oscillations or pulse waves. Figure S1 shows a schematic for a self-oscillating gel driven by pulse waves under periodic boundary conditions. In fig. S2, we examine the relative stability of different wave numbers in the absence of illumination by simulating 1000 independent runs with different random initial conditions at gel sizes L = 126.0, 315.0, and 630.0. At all gel lengths studied, the most probable outcome is a single wave propagating over the entire gel. To speed up the calculation, we initialize a pulse wave in the dark with K = 1 in the subsequent simulations. The illumination effect on the locomotion of the system is studied by alternately allowing the wave to reach a stable behavior, recording its characteristics, and then increasing the illumination by 0.001. For all spatiotemporal plots, the data spacing in time and x are 0.5 and 1.0 U, respectively.

We define the curvature κ(x) at each grid point in the gel to characterize the local distortion of the chemical wave, which can be written asκ(x)=d2tdx2(1+(dtdx)2)3(11)where (x, t) gives the position of a wave peak. To simplify the calculation, we calculate the local κ(x) across each wave and then store the data in one row, which is assigned to the time t at which the wave peak crosses the gel center.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We are grateful to the Advanced Analysis and Computation Center of CUMT for the award of CPU hours to accomplish the computations in this paper. Funding: This work was supported by the National Natural Science Foundation of China (grant nos. 21573282 and 21972165), the Natural Science Foundation of Jiangsu Province (grant nos. BK20171186), and the NSF (grant no. CHE-1856484). Author contributions: All authors contributed to discussion, the execution of the study, and theoretical analysis. L.R., Q.G., and I.R.E. conceived the research and wrote the paper. L.R. and L.Y. developed the numerical codes. L.R., L.Y., Q.G., R.T., and J.W. performed the calculation and prepared the figures and movies. Q.G. and I.R.E. supervised the research. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article