The Symmetry Basis of Pattern Formation in Reaction-Diffusion Networks
Ian Hunter, Michael M. Norton, Bolun Chen, Chris Simonetti, Maria Eleni Moustaka, Jonathan Touboul, Seth Fraden
TThe Symmetry Basis of Pattern Formation in Reaction-Diffusion Networks
Ian Hunter, ∗ Michael M. Norton, ∗ Bolun Chen,
3, 4
Chris Simonetti, Maria Eleni Moustaka, Jonathan Touboul,
3, 5 and Seth Fraden † Brandeis University Physics, Waltham MA, 02453 USA Center for Neural Engineering, Department of Engineering Science and Mechanics,The Pennsylvania State University, University Park, PA 16802 USA Brandeis University Volen National Center for Complex Systems, Waltham MA, 02453 USA Department of Physics Boston University, Boston MA, 02215 USA Brandeis University Mathematics Department (Dated: March 2, 2021)In networks of nonlinear oscillators, symmetries place hard constraints on the system that canbe exploited to predict universal dynamical features and steady-states, providing a rare genericorganizing principle for far-from-equilibrium systems. However, the robustness of this class of the-ories to symmetry-disrupting imperfections is untested. Here, we develop a model experimentalreaction-diffusion network of chemical oscillators to test applications of this theory in the context ofself-organizing systems relevant to biology and soft robotics. The network is a ring of 4 identical mi-croreactors containing the oscillatory Belousov-Zhabotinsky reaction coupled to nearest neighborsvia diffusion. Assuming perfect symmetry, theory predicts 4 categories of stable spatiotemporalphase-locked periodic states and 4 categories of invariant manifolds that guide and structure transi-tions between phase-locked states. In our experiments, we observed the predicted symmetry-derivedsynchronous clustered transients that occur when the dynamical trajectories coincide with invari-ant manifolds. However, we observe only 3 of the 4 phase-locked states that are predicted for theidealized homogeneous system. Quantitative agreement between experiment and numerical simu-lations is found by accounting for the small amount of experimentally determined heterogeneity.This work demonstrates that a surprising degree of the network’s dynamics are constrained by sym-metry in spite of the breakdown of the assumption of homogeneity and raises the question of whyheterogeneity destabilizes some symmetry predicted states, but not others.
I. INTRODUCTION
Network science unifies the study of disparate physicalsystems that can be cast as discrete sets of interactingdynamical units [1]. Here, we focus on networks of self-driven oscillators for which this simple framework pro-vides profound insights into systems ranging from elec-trical power grids to biological neural networks known ascentral pattern generators (CPG) responsible for coordi-nating autonomous animal locomotion [2–6].The design of networks that generate bespoke spa-tiotemporal patterns is a great challenge because univer-sal organizing principles for far-from-equilibrium systemsare exceedingly rare. Exploiting network symmetry isone way to meet this challenge. Symmetries place hardconstraints on the network dynamics of self-driven os-cillators by dictating that certain transient features andsteady-state patterns must exist. Specifically, the the-ory of equivariant dynamics describes how the symme-tries of the network affect the symmetry of the networkdynamics[7]. For a given symmetry operation, such asa permutation of network nodes, some sets of points instate-space will be unchanged and consequently, their dy-namics must also remain the same [7, 8]. ∗ These two authors contributed equally † To whom correspondence should be ad-dressed:[email protected]
A class of results derived from group theory arises bycombining the spatial symmetry of the network with thetemporal symmetry of the oscillators, providing a naturalframework for describing spatiotemporal patterns[3, 7–11]. One, the H/K theorem, allows enumeration of allsymmetry derived patterns in which phase-locked nodesco-evolve because they receive the same input from theirneighbors[9]. Remarkably, some of the predicted patternsare far from obvious and bear little resemblance to thegeometric symmetry of the network. Significantly, thesepatterns are universal. They depend only on the cou-pling topology and are independent of all system specificdetails regarding the nature of the non-linear oscillatorsthemselves and even whether or not the coupling is non-linear. However, these striking results derive from thestrong assumption that classes of nodes in the networkand their interconnections are strictly identical [8, 11]Golubitsky and colleagues applied the H/K theorem,along with a few plausible assumptions, to make asurprising prediction in neuroscience; the minimal net-work architecture of central pattern generators in allquadrupeds can be determined simply by cataloguingthe aggregate of gaits observed across species. Or,in other words, that form follows function in neuronalnetworks[3]. The existence of such CPGs is controversialin the case of mammals, but evidence exists for otherorganisms [12–14].In this work, we experimentally study oscillatory chem-ical reaction-diffusion networks and examine the dy-namics through the lens of symmetry-based network a r X i v : . [ n li n . PS ] F e b theories[8, 15]. The significance of studying a self-contained reaction-diffusion system lies in the potentialfor fabrication of autonomous devices that organize theirspatiotemporal dynamics through processes analogous toliving systems.Our goal here is to ascertain whether symmetry canserve as a conceptual basis and engineering principle forthe structuring of spatiotemporal patterns on chemicalnetworks, which can be equally applied to understand-ing biological neural networks and engineering chemicalnetworks for soft robotics[16].The title of this paper is a paean to Turing who wasthe first to consider the theory of pattern formation indiscrete reaction-diffusion networks [17]. Although themajority of Turing’s paper, “The chemical basis of mor-phogenesis,” focused on static, spatially varying patterns,Turing also predicted spatiotemporal pattern formation,including standing and traveling chemical waves, whichhave been observed in chemical networks[18]. It was thislatter aspect of Turing’s theory that motivated us to ex-perimentally test symmetry-based network theory usingcoupled chemical oscillators. However, because Turing’slinear stability analysis is limited to the onset of patternformation, we were motivated to employ the symmetryapproach of the H/K theory because it is universal, hold-ing true for non-linear oscillators and for all time.To test the network theory, we develop a minimallycomplex experimental system consisting of a ring of 4identical, nanoliter sized, chemical reactors containingthe Belousov-Zhabotinsky (BZ) oscillating reaction andcoupled to nearest neighbors by diffusion. This experi-mental system oscillates stably for about 70 periods[19],which is an order of magnitude longer than reported inprevious studies of self-organized networks[4, 20, 21]. Tothoroughly explore state-space, we perform hundreds oftrials by running experiments simultaneously on multiplecopies of the network resulting in an order of magnitudegreater number of experiments than done previously withdifferent chemical networks[16]. We model this reaction-diffusion network at two levels of description. The mostdetailed is a mathematical model of the system explicitlydescribing the BZ reaction chemistry, which we assumeoccurs only in the reactors, with coupling between near-est neighbors caused by diffusion of a subset of the BZchemicals through the intervening PDMS. We theoreti-cally reduce this reaction-diffusion network into a simplerphase model to analyze the predictions of the theory ofequivariant dynamics. These models contain far fewerfree parameters than independent measurements. This,the large ensemble of experiments and their longevity al-lows quantitative comparisons between theory and exper-iment leading to firm conclusions regarding the applica-bility of idealized network dynamics to model the steady-state and transient dynamics of this self-organized sys-tem in which both the oscillators and coupling are fullychemical.Our experiments reveal an intricate array of tran-sient and phase-locked spatiotemporal chemical dynam- xy BZ Inhibited BZ PDMSActinic Light + BZ PDMS1 234
Time/T o (a) (b)(c) (d) Theory
RFRHLHLF I n t en s i t y [ A . U ] FIG. 1. (a)
Schematic of a network of a ring of 4 inhibitorycoupled oscillators. Indexing of nodes is indicated as eithera number (1,2,3,4), or leg of a quadruped (LF, RF, RH, LH)with L left, R right, F front and H hind. (b)
Schematicof the experimental system. The reactors are divots in thePDMS, filled with BZ and sealed between 2 glass plates. (c)
Photograph of BZ filled 4-ring network. Actinic light illumi-nates BZ in a channel surrounding the network and providesa constant chemical boundary condition. (d)
Two adjacentreactors (red and blue traces) in the network oscillating 180 ◦ out-of-phase with each other. Top:
Measured transmittedintensity versus time.
Bottom:
Simulated oxidized catalystconcentration [mM] versus time. In both, time is rescaled byoscillation period T , indicated by the two arrows. ics. However, when we reduce these complex non-stationary solutions into the state-space of phase re-lationships, the H/K theorem allows us to representhigh-dimensional chemical dynamics in terms of sim-ple, model-independent geometric objects in the formof planes, lines, and points that are readily visualized.This geometric perspective leads to an appreciation ofaspects of the experimental behaviors that are directlyimposed by symmetries, thereby providing conceptualunderstanding of the complex dynamics, which comple-ments and enriches the quantitative comparison betweentheory and experiment. II. RESULTSA. Experimental Reaction-Diffusion Network
We designed a reaction-diffusion network consistingof a ring of four diffusively coupled nanoliter volumebatch reactors laid out in a square 2x2 lattice with near-est neighbor coupling (see Fig. 1). Previously, we em-ployed emulsions containing the BZ oscillating reactionto study reaction-diffusion networks[18, 22–28]. But thediffusive coupling between surfactant stabilized emulsiondrops was difficult to characterize and manufacturing ofthe networks was challenging, which contributed to alarge degree of variability between experiments. Here,to improve reproducibility we manufactured these reac-tors to high precision from elastomeric PDMS using softlithography techniques and filled the reactors with theoscillatory BZ reaction as described previously[16, 19],illustrated in Fig. 1 and in Appendix II C. To obtaina large statistical sample of trajectories we made de-vices that combined 9 or 16 copies of the 2x2 network[16]. To optimize homogeneity in the chemical concen-trations of each of the reactors, we simultaneously filledthe entire set of networks by pipetting a drop of BZ thatfloods all the reactors before sealing the sets of reac-tors by clamping the PDMS between two glass plates[Supplementary Material [29] Fig. 1-2, SI videos in [16]].The chemical coupling between adjacent reactors arisesfrom the permeation of chemical species through the in-tervening PDMS wall and mainly consists of bromine-induced inhibition, with a weaker activator coupling, per-haps by bromous acid and the bromine dioxide radi-cal [16, 18, 22–24, 26–28, 30–34]. After mixing the BZreagents, pipetting them onto the PDMS networks, seal-ing the networks, and placing the sample in the dark foran induction period of 20 minutes, it was observed thatall reactors began to oscillate and collectively form spa-tiotemporal patterns [Fig. 1(d)] [16].The reactors form a closed system and consequentlythe oscillators have a finite lifetime as the reactantsare consumed and waste products accumulate. How-ever, although the amplitude of the chemical oscil-lations decreases over time, the oscillators maintaina nearly constant period for a duration of order 70oscillations[19]. Based on this long term stability, weassume that the underlying phase dynamics of the indi-vidual BZ oscillators remains constant during the dura-tion of the experiment, thus allowing us to study phaserelationships between reactors as they evolve over time[Supplementary Material [29] Fig. 4, Movies S1-4]. Each4-ring network is isolated from the environment becausethe reactors are surrounded by a zone of photosensitiveBZ that is held at constant chemical conditions by theapplication of actinic light[16]. We also assume that eachchemical reactor is well mixed, ignoring any spatial vari-ation of chemical concentrations, because the size of thereactor is small compared with the length scale of diffu-sion, e.g. w < √ Dτ with w the width of each squarereactor ( w = 62 µ m), D , the diffusion constant of eachBZ chemical ( D ∼ − m s − ) and τ , the duration of aBZ oscillation ( τ ∼ B. Theory and the Role of Symmetry
The fullest description of the dynamics of the 4-ringnetwork that we consider is a reaction-diffusion networkmodel. It focuses on the time dependent concentrationsof the well mixed chemicals in each reactor, denoted as(¯ c ( t ) , ¯ c ( t ) , ¯ c ( t ) , ¯ c ( t )), where ¯ c j ( t ) is a vector of concen-trations in the j th reactor with indices as in Fig. 1A. As-suming that the reactors are identical, the behaviors areexpected to have the same symmetries as a square, e.g. 3 rotations and 4 reflections. Associated with this symme-try group, the H/K theorem predicts invariant manifolds,subspaces in which the dynamics remains confined, thatare universal to any ring of 4 oscillators, enumerated inTable I. Although the theory is more general, we restrictourselves to the case in which all the nodes are on thesame limit cycle, as this corresponds to experiment. Withthis assumption the H/K theorem guarantees any systemof 4 oscillators with square symmetry possesses 8 cate-gories of invariant manifolds, including 4 categories ofphase-locked periodic states. These states are thereforeefficiently described by considering the phase relationshipbetween pairs of reactors, defined as the fraction of pe-riod they are shifted from each other on their commonlimit cycle.Invariant manifolds are denoted by a pair of symme-try operations, (H,K). The first symmetry, H, representsan exchange of nodes that results in the same state sub-ject to one or more phase-shifts. The second symmetry,K, indicates symmetries under which the system is un-changed.Depending on the constraints imposed by H and K,these solutions may either maintain fixed phase relation-ships among all four nodes resulting in phase-locked so-lutions, or leave 1- or 2-dimensional freedom on these re-lationships, as enumerated in Table I. The 4 categories ofphase-locked periodic states are spatiotemporal periodicpatterns that correspond to 6 point invariant manifoldsin the phase difference space that can be identified withgaits of quadrupeds enumerated in Table I and visualizedin Fig. 2(b). The first two categories are Pronk in whichall the legs advance simultaneously and
Trot for whichdiagonal legs are in phase, and the two diagonal pairsof legs are half a period out of phase.
Pace and
Bound form one category and we refer to them interchangeablyin the remainder of the text. In Pace, legs on each sideare in phase and opposite sides out-of-phase, while forBound, legs on opposite sides are in phase and the frontlegs out-of-phase with the hind legs.
Clockwise (counterclockwise) Rotary Gallop is another category in whichthe legs advance in a clockwise (counter clockwise) man-ner with each leg advancing a quarter of a period laterthan the preceding leg.The remaining 4 categories correspond to higher di-mensional invariant manifolds (lines or planes) that con-tain trajectories maintaining partial symmetries. Along1-dimensional linear invariant manifolds, the network canbe split into two pairs of reactors, such that within pairsthe reactors are in phase or antiphase, while betweenpairs reactors have an arbitrary phase-shift. Along 2-dimensional manifolds, two nodes oscillate in phase andthe other two nodes are at arbitrary phase-shifts. In fact,the 2-dimensional manifolds intersect the 1-dimensionalmanifolds, and the 1D manifolds intersect the phase-locked 0-dimensional manifolds [Fig. 2(a)]. Heuristically,these higher dimensional manifolds often act as privilegedpathways that both guide and structure transient transi-tions between the phase-locked states. Beyond predictingthe existence of these invariants, the H/K theorem nei-ther prescribes their stability nor precludes the existenceof others. To address questions of stability and existenceof additional manifolds requires a specific model of theoscillators and their connections.
C. Chemical Kinetic and Phase Models
To model the reaction-diffusion dynamics of our exper-iments we use the Vanag-Epstein model of the BZ reac-tion, which treats the chemical kinetics of 4 BZ chemicals,Br – , HBrO , Ferroin, and Br , combined with a diffu-sive coupling term between chemical species in which thecoupling strength is fitted from the data. Noting that¯ c i ∈ R are the concentrations of these 4 chemicals inreactor i ∈ { , , , } and R : R R denotes theVanag-Epstein BZ model vector field [16, 18, 22–24, 26–28, 30–34], we obtain the equation: ddt ¯ c i = ¯ R (¯ c i ) + X j =1 A ij µ (¯ c j − ¯ c i ) (1)where µ ∈ R × accounts for the chemical coupling ma-trix between two adjacent cells, and depends on the per-meability of the PDMS for each chemical species and thegeometry of the reactors, while A ∈ R × denotes the ad-jacency matrix between two reactors that is determinedby the network topology [Appendix II H for details onthe model, choice of parameters and fits of free parame-ters]. This model ignores spatial concentration gradientsinside the reactors, effectively treating reactors as pointscorresponding to nodes of the network. The model alsoneglects the occurrence of chemical reactions within thePDMS, which acts as a connector that couples adjacentnodes.Under the assumption that all reactors are oscillatingon the same limit cycle, we can parameterize the time de-pendent concentrations through the phase of that cycle,as proposed by Winfree [35] and widely used in variousapplications [28, 36–40]. In this abstraction, the phasevariable naturally progresses linearly from 0 to 2 π at afrequency ω = 2 π/T . Perturbing the chemical concen-tration of an individual reactor will lead to a modificationof the phase that depends on the chemical species thatis perturbed and phase of the reactor; this function iscalled the phase response curve (PRC). The impact onthe phase of one reactor due to diffusive coupling froma neighbor can be summarized through the interactionfunction, H , derived from convolving the PRC with thediffusive coupling between reactors and averaging over aperiod [37] [Appendix II H]. Notably, this reduction ofthe chemical model of Eq. 1 to a phase model introducesno new parameters. Best fits between experiment and TABLE I. Symmetry required invariant manifolds for an os-cillator network possessing square or Dihedral 4 ( D ) sym-metry. D , all symmetries of a square; D pn , reflection across n diagonals; D sn , reflection across n , vertical or horizontal,axes; Z , 90 ◦ rotation; Z , 180 ◦ rotation; 1, no operation.The first 4 classes of manifolds are phase-locked states. Thecolumn marked “Phase” graphically indicates the spatiotem-poral pattern with symbols representing the phase in percent-age of the period T , white circle - 0%; white/black - 25%;black circle - 50%; black/white - 75%. T denotes the periodof each oscillator in a given invariant manifold and can varyfrom manifold to manifold. The second 4 classes of mani-folds are symmetrically clustered states, related by arbitraryphase shifts f , f , which vary from 0 to 1 as fraction of aperiod. The graphical representation of nodes in the column“Phase” have solid, striped, or dot motifs. Different motifsare related by an arbitrary phase shift. Similar motifs withopposite background colors are antiphase with each other.Point invariant manifolds: Name(H,K) Phase ¯ c ¯ c ¯ c ¯ c Pronk ( D ,D ) ¯ c ( t ) ¯ c ( t ) ¯ c ( t ) ¯ c ( t ) Trot ( D ,D p ) ¯ c ( t ) ¯ c ( t + T ) ¯ c ( t ) ¯ c ( t + T ) Pace ( D s ,D s ) A Bound ( D s ,D s ) B ¯ c ( t )¯ c ( t ) ¯ c ( t + T )¯ c ( t ) ¯ c ( t + T )¯ c ( t + T ) ¯ c ( t )¯ c ( t + T ) CW Gallop ( Z , A CCW Gallop ( Z , B ¯ c ( t )¯ c ( t ) ¯ c ( t + T )¯ c ( t − T ) ¯ c ( t + T )¯ c ( t + T ) ¯ c ( t − T )¯ c ( t + T )Linear invariant manifolds: ( D s ,D s ) A ( D s ,D s ) B ¯ c ( t )¯ c ( t ) ¯ c ( t + f T )¯ c ( t ) ¯ c ( t + f T )¯ c ( t + f T ) ¯ c ( t )¯ c ( t + f T ) ( D s , A ( D s , B ¯ c ( t )¯ c ( t ) ¯ c ( t + T )¯ c ( t + ( f + ) T ) ¯ c ( t + f T )¯ c ( t + f T ) ¯ c ( t + ( f + ) T )¯ c ( t + T ) ( Z , ¯ c ( t ) ¯ c ( t + f T ) ¯ c ( t + T ) ¯ c ( t + ( f + ) T )Planar invariant manifolds: ( D p ,D p ) A ( D p ,D p ) B ¯ c ( t )¯ c ( t ) ¯ c ( t + f T )¯ c ( t + f T ) ¯ c ( t )¯ c ( t + f T ) ¯ c ( t + f T )¯ c ( t + f T ) model are obtained with the interaction function H thatarises from a combination of Br and HBrO , as shownin Fig. 6(a). This leads to the phase equation: ddt φ i = ω + k X j =1 A ij H ( φ j − φ i ) (2)where k the diffusive coupling rate.Noting that the right hand side of Eq.2 depends onlyon phase difference θ ij ≡ φ i − φ j ; we therefore arbitrar-ily choose the three phase differences ¯ θ = ( θ , θ , θ )as the new system variables and recast the dynamics ac-cordingly, ddt ¯ θ = ¯Ψ(¯ θ ) (3)with ¯Ψ(¯ θ ) following directly from Eq. 2. As each of thephase differences is periodic on (0 , π ], the state-spaceis a 3-torus. Although the 3-torus cannot be drawn inthree dimensions, it is equivalent to a Cartesian cubewith periodic boundaries, allowing visualization of thefull dynamics. Dynamics in the State-space of Phase Differences
This new coordinate system transforms the invariantmanifolds identified by the H/K theorem in Table I tosimple, geometric objects: points, lines, and planes, enu-merated in Table II. This transformation enables theconsequences of the H/K theorem on the dynamics instate-space to be visualized in a way that would be im-possible in the full chemical model, Eq. 1. The pointinvariant manifolds Pronk, Bound, Trot, Rotary Gallopbecome steady-states in this new frame, rather than high-dimensional limit cycles. We are able to readily classifythem as either attractors, repellers or saddles, accordingto whether the velocity vectors surrounding the steady-state point inward, outward or change sign depending onorientation, respectively [41]. Beyond the steady-statesin the form of points, Fig. 2A shows a state-space struc-tured by additional invariant manifolds in the form oflines and planes, all required by the spatiotemporal sym-metries of the oscillator network.The dynamical system of Equation [3] predicts thatthe network is multistable, with the point H/K mani-folds forming competing attractors. We simulated ex-haustively the model and observed that each initial condi-tion flows to one of the four categories of phase-locked at-tractor states; (1) Pronk, (2) Pace/Bound, (3) Trot, and(4) CW/CCW Gallop [Fig. 3(a)(b)]. Furthermore, weshow in the Supplementary Material [29] Sec. IIB thatthe system predicts the six phase-locked states are lin-early stable, and thus attractors.We found that the theoretical model also possessesmany unstable, saddle phase-locked steady-states in ad-dition to the six H/K derived attractors. In fact, topol-ogy predicts the existence of unstable steady-states.The topological index of both attractors and saddleswith one attracting direction is +1, while the indexof both repellors and saddles with two attracting di-rections is -1. Topology requires that the sum of thetopological indices of all the steady-states must equal0, the Euler characteristic of a 3-torus, as shown inthe Supplementary Material [29] Sec. IV. Given the sixH/K point invariant manifolds are attractors, we con-clude there must be at least six unstable steady-stateslocated in the 3-torus to satisfy the required charge neu-trality. Moreover, unstable states organize the separa- trices between invariant manifolds containing more thanone attractor. We numerically searched for the requiredunstable states and found 158 saddle and 4 unstablesteady-states dispersed throughout the state-space shownin Supplementary Material [29] Fig. 7. A notable andunexplained fact is that of the 168 numerically identifiedphase-locked steady-states, the sole attractors are the sixinvariant point manifolds required by the H/K theorem.
TABLE II. Symmetry required invariant manifolds parame-terized by relative phase. f , f , which vary from 0 to 2 π . Allrepresentations, modulo 2 π , are shown.Point invariant manifolds: Name Phases Hyperplane
Pronk ( θ = θ = θ = 0) Trot ( θ = θ = θ = π ) PaceBound ( θ = − θ = π, θ = 0)( θ = − θ = 0 , θ = π ) CW GallopCCW Gallop ( θ = θ = θ = + π )( θ = θ = θ = − π )Linear invariant manifolds: ( D s ,D s ) A ( D s ,D s ) B ( θ = − θ = f , θ = 0)( θ = θ = 0 , θ = f ) ( D s , A ( D s , B ( θ = π, θ = f , θ = − θ )( θ = f , θ = π, θ = π ) ( Z , ( θ = f , θ = π − f , θ = f )Planar invariant manifolds: ( D p ,D p ) A ( D p ,D p ) B ( θ = − θ , θ = f , θ = f )( θ = f , θ = f , θ = − θ ) We numerically determined the basins of attraction ofeach attractor by dividing the 3-torus into a fine gridand identifying each initial point with the attractor towhich it flowed, as shown in Fig. 3A-B [42]. The Pronk,Bound and Rotary Gallop basins are smooth, closed vol-umes while the Trot basin fills the rest of the state-space[Fig. 3(b)]. The state with largest basin of attractionis Trot, followed by Rotary Gallop, Bound and Pronk.The attraction basin of the Bound state is anisotropicand aligned with the ( D s , D s ) invariant manifolds [Fig.3(a),4(a)]. Fig. 3(a) and 4(b) reveal that trajectories re-main near the ( D p , D p ) invariant manifolds as they flowto Trot. Theory predicts that the network’s trajectoriesflowing towards its attractors are constrained and shapedby H/K linear and planar invariant manifolds.To further elucidate how symmetric invariant man-ifolds guide and structure dynamics, we focus on thetransverse dynamics, namely flows perpendicular to the (D p1 ,D p1 ) B (D s1 ,D s1 ) A (D s1 ,1) A (Z ,1)(D s1 ,D s1 ) B (D p1 ,D p1 ) A (D p1 ,D p1 ) B (D s1 ,1) B Pronk TrotBoundCWRotary Gallop CCWRotary Gallop(D p1 ,D p1 ) A Pace (a) Theory: Invariant manifolds
TrotPronk
Pace
Other S pa t i o t e m po r a l pa tt e r n s (.45,0.,.46) (.5,.5,.5)(0.,0.,0.) (.5,.35,.45) O b s e r v ed
26 11 73 76 (c) Experiments: Steady-states ( , , ) (b) Theory: Stable steady-states C l u s t e r s S pa t i o t e m po r a l pa tt e r n (0,0,0) T i m e ( , , ) RotaryGallopTrotPronk
Pace
FIG. 2. (a)
The 8 categories of invariant manifolds for a network of 4 nodes with square symmetry predicted by the H/Ktheory are presented in the state-space of phase differences. There are 4 categories of phase-locked periodic states (points), 3categories of lines and 1 category of planes. The state-space is periodic and invariant manifolds that are identical mod 2 π arerendered in translucent colors. (b) (First row) The 4 point manifolds predicted by theory, Pronk, Pace, Trot and Gallop, arerepresented as space-time plots, (second row) as networks with fixed phase differences and (third row) as a triplet of phasedifferences, ( θ , θ , θ ), with units of fraction of a period. (c) (First row) In experiment, 3 of the 4 phase-locked gaits areobserved, while Rotary Gallop is never observed. However, a myriad of states with an unclear correspondence to theory areobserved, labeled “Other”. The “Other” states are more similar to Trot than Pronk, Pace or Gallop in terms of space-timeplots and phase differences. However, Trot and “Other” are qualitatively distinct states because in Trot diagonal nodes oscillatein-phase while in “Other” they oscillate with large phase shifts between them, shown by the red and green bars over the space-time plots. (Second row) The number of times a state is observed is recorded and (Third row) the measured phase differenceis listed. Videos of experiments shown in Movies S1-S4. invariant manifolds. We combined the unique systemdynamics, Eqns. 2 and 3, and the universal struc-ture of the line and plane invariant manifolds to semi-analytically compute the manifolds’ transverse Lyapunovexponents, which measure the local attraction or re-pulsion rate of trajectories to or from them [43–45],in terms of the interaction function, H [AppendixII L]. The majority of the domains of the ( D s , D s ) and( D p , D p ) manifolds are covered with negative exponents[Fig. 4(c)]. This causes trajectories to collapse and re-main on these invariant manifolds [Fig. 4(a)(b)] [45].While the regions with negative exponents organize theflows along the invariant manifolds, the small regions ofpositive exponents contain the separatrices of attractorson the same invariant manifolds, shown as saddles in theSupplementary Material [29] Fig. 7. In this way, the the-ory combines the restrictions of symmetry and the uniquesystem dynamics to predict both the basins of attractionof the attractors and the symmetric, clustered transienttransitions along the linear and planar invariant mani-folds that connect the attractors.The higher order H/K invariant manifolds are notthe only invariant manifolds controlling transient flows.Additional invariant manifolds arise in the vicinity ofsteady-states because linear stability analysis permits aneigenvector decomposition of the dynamics. We see such an invariant manifold exists about the Rotary Gallop at-tractor as all the surrounding trajectories coalesce into aplane [Fig. 3(a)]. This invariant manifold is not alignedwith any H/K invariant manifold and arises from theunique system dynamics rather than the H/K theoremor symmetry. Experimental Observations of Dynamics
To compare the geometric framework provided by thetheory of equivariant dynamics with experiment, we de-fined a reactor’s phase in reference to the moments ofmaximum oxidation, corresponding to a maximum intransmitted light intensity [Fig. 1(d)]. These momentswere attributed phase 2 π , and phase was defined as thefraction of time spent between two peaks. This allows usto directly compare theory with experiment by measuringthe experimental phase differences between 3 sequential,adjacent pairs of reactors ( θ , θ , θ ), and to charac-terize the dynamics relative to invariant manifolds andtheir stability derived from the model.We conducted 318 experiments of which 186, or 58%,phase-locked as defined in Appendix II H. During theinterval of time soon after the first oscillation and be-fore phase-locking, the oscillation periods of all of the (a) Simulated trajectories Rotary Gallop TrotBound PacePronk ExperimentSimulationPronkBasin BoundBasin Rotary GallopBasin TrotBasinPaceBasin TrotBasinPronk TrotBoundOther BoundBasinRotary GallopBasin N: (c) Basins of attraction: Experiment vs theory(b) Basins of attraction: Experiment vs theory FIG. 3. Basins of attraction. States are labeled as in Fig. 2. (a)
Simulations of eqn 3 show that all trajectories converge tothe H/K point invariant manifolds, Pronk, Bound, Rotary Gallop or Trot, depending on initial condition. Video of differentperspectives in 3D are shown in movie S5. The corresponding plot of all experiment is shown in movie S6. (b)
Basin ofattraction for experiments and theory. The theoretical basins of attraction are colormaps computed from 5399 simulations,including those in (a). For example, if an initial state is colored red, then it will flow to the Trot attractor. The experimentalbasins of attraction are partially reconstructed by disks located at their initial condition and colored by the attractor to whichthey converge, with the color code indicated in (c). We artificially shift slightly the experimental points with nearby initialconditions so that they are revealed. Video of different perspectives in 3D are shown in movie S7. (c)
The probability, P ( A F | A I ) , that an experiment which started at a set of phase differences in a given theoretical basin, A I , converges to eachof the 4 experimental attractors, A F . The number of observations is listed in the row labeled “N”. For example, of the 2 initialconditions corresponding to Rotary Gallop, one ended up in the Trot basin and one in the Other basin, but neither went tothe predicted basin. In contrast, of the 148 states initially in Trot, 48.7% went to Other, 2.0% to Pronk, 1.3% to Bound and48.0% to Trot. reactors remained similar, ± Other if theyare more than 1.0[rad] from each of the H/K attractors.Of the observed steady-states, a majority of 59% cor-respond to the predicted Pronk, Bound and Trot pointinvariant manifolds [Fig. 2(c)]. We were initially baffledby the remaining 41% of the observed phase-locked statesas they are located a distance from each of the 4 classesof attractors that exceeds the aforementioned thresholdradius [Fig. 2(c)], raising the question of their origin asthey were not predicted by theory. Strikingly, Fig. 4(a)(b) shows that experimental tra-jectories starting near invariant manifolds ( D s , D s ) and( D p , D p ) closely follow the dynamics predicted by thetheory. Such a consistency between experimental obser-vations and clustered transient states predicted on thebasis of network symmetry alone is quite remarkable.This suggests that the universal properties dictated solelyby the symmetry of the system not only predict station-ary behaviors of an experimental system, but also con-strain transient dynamics from an initial condition to astationary state.There are three other noteworthy comparisons to makebetween theory and experiment. Firstly, we find experi-ments and theory have similar shaped Pronk, Bound andTrot basins of attraction, as illustrated in Fig. 3(b)(c).Secondly, the theoretically predicted Trot attractor issymmetrically surrounded by an extended cloud of phase-locked states, denoted Other, which are significantly farfrom Trot, yet within the predicted Trot basin of attrac-tion. These observations suggest that the experimen-tal Other states are associated with the predicted Trotstate, as shown in Fig. 3(b)(c) and Fig. 5(a). Thirdly,the Rotary Gallop states were absent in all our experi-ments. This is particularly surprising because the theoryattributes to that state a large basin of attraction. =- In-plane dynamics2000Time[s] 3000 5000SimulationExperimentPronkPronk TransverseLyapunov exponent [1/ks]-1.8-.9103.46.7TheoryExperimental TheoryExperimentalTransverseLyapunov exponent [1/ks]-1.8-.9103.46.7Time[s] 1000 2000 3000SimulationExperimentTrotPacePace Pronk TrotTrotPronk= In-plane dynamics TransverseLyapunov exponent [1/ks]-1.8-.9103.46.7 (a)Dynamics on1D invariant manifold (b)Dynamics on 2D invariant manifold (c)Theory: Dynamics on invariant manifolds Pronk(D s1 ,D s1 ) A (D p1 ,D p1 ) A Pace Trot (A) (D s1 ,D s1 ) A Pronk(D p1 ,D p1 ) A (B) PronkPronk Pace
FIG. 4. Transient dynamics along higher order H/K invariant manifolds in experiment and simulation. (a)
Space-timeplots from an experiment and simulation with a near-Pronk initial state transitioning to Pace. States form 2 symmetricclusters corresponding to the ( D s , D s ) invariant manifold. In the lower panel the experimental trajectory is shown as an arrowtraveling through a 2D slice of state-space superimposed over the theoretical velocity field. Video of experiment synchronizedto progression along space-time plot and trajectory in state-space shown in movie S8. (b) Space-time plots from an experimentand simulation with a near-Pronk initial state transitioning to Trot. The transition corresponds to the ( D p , D p ) invariantmanifold. Video of experiment synchronized to progression along space-time plot and trajectory in state-space shown in movieS9. (c) The invariant manifold surfaces attract or repulse in a state-dependent manner. The analytically computed transverseLyapunov exponents [Appendix II L] are shown via heatmaps. When positive it indicates nearby trajectories are repulsed fromthe invariant manifold. When negative it indicates attraction. Both the 2D ( D p , D p ) and 1D ( D s , D s ) invariant manifolds arelargely attracting. Video of 3D perspective of plot in movie S10. Inclusion of Slight Heterogeneity in Theory
A hypothesis to account for these discrepancies be-tween theory and experiment is to consider hetero-geneities. Indeed, chemical systems differ from theoreti-cal models in that they are bound to be heterogeneous,and in particular imperfectly symmetrical [1, 37]. For in-stance, modeling BZ micro-oscillators with small degreesof heterogeneity in reactor chemistry (in turn associatedwith heterogeneous frequencies), or in reactor volume,led in various situations to better match experimental re-sults [18, 26, 28]. We thus tested whether small degreesof heterogeneity between reactors could indeed explainpart of the discrepancy between theory and experiments.To this purpose, we defined a variation of our networkmodel including heterogeneities between reactors, which,using phase reduction [see detail in Appendix II H] led usto analyze a phase model with heterogeneous frequencyof type[38]: ddt φ i = ω i + k X j =1 A ij H ( φ j − φ i ) (4)To fit equation 4 to a trajectory of phase differ-ences measured during an experiment requires a spe-cific set of non-zero unperturbed frequency differences ω i − ω j [Appendix II H]. Although the best fit un-perturbed frequencies are different for each experiment,their statistical distribution fits a Laplacian probabil-ity distribution function corresponding to unperturbedfrequencies ω i having a percent coefficient of varia-tion of ± . N = 5 , . N =34 , N = 186)in experiments. This 23-fold decrease in percentage ofRotary Gallop steady-states is significantly larger thanthe 1.5, 2.3 and 3.2-fold decreases for Trot, Bound andPronk. Further, the Trot, Bound, and Pronk states insimulation are clustered in a manner corresponding toexperiment [Fig. 5(b)]. In particular, we observe states,termed Other, that form a large cluster of phase-lockedstates centered about Trot [Fig. 5]. We therefore con- (a) Experimental steady-states (b)
Heterogeneous simulationsteady-states
PronkBound TrotOther Pace RightRotary Gallop LeftRotary GallopPronkBound TrotOther Pace
FIG. 5. Comparison between experiment and simulationwith heterogeneous oscillator frequencies. Phase-locked statesare colored by which H/K point manifold they are closest to,or as “Other” if they are more than 1.0 [rad] away from eachgait. States are labeled as in Fig. 2. (a)
Experiments. (b)
Simulations with heterogeneity. A 3% percent coefficient ofvariation in unperturbed frequencies, ω i , caused a 23-fold re-duction in observed Rotary Gallop steady-states compared tothe homogeneous simulations in Fig. 3(a)(b). Videos of 3Dperspectives of plots shown in movies S11 and S12, respec-tively. clude that small levels of heterogeneity can indeed ac-count for the main discrepancies between the theory andthe models. DISCUSSION
Symmetry principles have been used to design net-work topologies of electronic oscillators that generate de-sired dynamics [5, 6]. The correspondence between the-ory and experiment was ascribed to the precision of thefabrication process of these electromechanical devices,which possessed intrinsic frequencies with only 0.001%percent variation, ensuring the equivalence of each net-work node and connection [6]. While these experimentsdemonstrated the relevance of equivariant dynamics toexperiments, it is unresolved whether the predictions ofH/K theorem persist in the biological milieu which neverachieve such interchangeability [35].Approaching biological systems directly in a similarmanner is difficult because the intrinsic dynamics areboth complex and unknown. One study found that os-cillatory slime mold confined to networks with specifiedsymmetries support aspects of the H/K theory [4]. Thebiological oscillators in this system possessed intrinsic fre-quencies with 10% percent variation [46]. Despite being4 orders of magnitude more heterogeneous than the elec-tronic system, symmetric dynamics emerged. Still, whilethese results are intriguing, because the underlying chem-ical dynamics of slime mold are not fully understood andthe total number of oscillations are small, it is impossi-ble to assess the stability of observed spatiotemporal pat- terns, transient dynamics, or the impact of heterogene-ity. Thus, minimally complex reaction-diffusion basedsystems provide an essential linkage between idealizedtheory and biology.Our experiments showed both consistencies and dis-crepancies with the symmetry-based theory. In particu-lar, most phase-locked states were recovered, and track-ing trajectories from an initial condition to a phase-locked state showed that the theory not only predictedthe steady-states, but also their basins of attraction and,more surprisingly, transient dynamics and transverse sta-bility along invariant manifolds.However, through combining theory and experimenton a model system, our results suggest that heterogene-ity eliminates some states, but not others, raising thequestion of understanding which states are the most sen-sitive to heterogeneity. A natural hypothesis would bethat the states disappearing upon addition of hetero-geneity would correspond to states which in the absenceof heterogeneity, would have weaker stability. Estab-lished theory concerning heterogeneity in oscillator net-works predicts the impact of heterogeneity is inverselyproportional to the system’s linear stability [28, 47].Following this approach, we computed the maximumLyapunov exponent for each of the phase-locked states.We found that Pronk and Bound have the same max-imum Lyapunov exponents, − × − [s − ], while Ro-tary Gallop has a threefold larger value, − × − [s − ][Supplementary Material [29] Fig. 6]. In other words, inthe absence of heterogeneity, the Rotary Gallop state hasa higher magnitude vector field pointing in towards itthan the other states. Therefore, interestingly, RotaryGallop has a stronger stability and a larger basin of at-traction relative to Pronk and Bound, and yet appearsmore sensitive to heterogeneity. This raises the questionof characterizing how symmetric states are affected byperturbations, a largely open theoretical question withimportant applications. CONCLUSION
Understanding how network structure controls spa-tiotemporal pattern formation remains a central prob-lem in network science. Analysis of spatial network sym-metries has led to great progress by illuminating mech-anisms behind the emergence of clustered, dynamicalstates. Specifically, tools for identifying group orbits [45]and equitable partitions [48, 49] have been particularlyfruitful in systematizing the identification of topology-required clustered states. Here, the low-dimensionalityof the representation of a 4-ring network by phase dif-ferences allows us to concretely illustrate complex butuniversal features of the dynamical landscape underpin-ning the emergence of clusters.We showed that a longstanding theoretical conjec-ture, that symmetry can dictate function in biolog-ical systems, can be used to rationally engineer a0spontaneously organizing reaction-diffusion network ofBelousov-Zhabotinsky oscillators. Thus, our results offerpromise for applying symmetry principles as a tool fordesigning out-of-equilibrium materials and understand-ing biological dynamics. Contemporary theory [3, 15]proves that oscillator networks with the same symmetryas the 4-ring we studied are required to share a universallist of phase-locked states and transient dynamics.By mapping both the high-dimensional chemical modeland experimental observations to a 3D state-space ofphase differences, we showed that the symmetry requiredinvariant manifolds take on simple forms. Point man-ifolds are phase-locked states with spatiotemporal pat-terns that we recognize as quadruped gaits. The higherdimensional manifolds consist of lines and planes in state-space that guide the transient dynamics from one pointmanifold to another. Additionally, network sub-clustersare sequentially synchronized when convergence to a fi-nal spatiotemporal pattern occurs along these manifolds.Thus the H/K theorem imposes a great deal of structureon the phase-locked and transient dynamics of the sys-tem that is dependent only on the network’s topologyand independent of any of the specifics of the oscillatorsand their coupling. These results therefore provide strongsupport to the hypothesis that symmetries in chemical orbiological neural networks organize functional patterns,such as locomotion [7, 15].An important aspect of this study is its exhaustiveness:analyzing a small network allowed a complete applica-tion of the H/K theorem. In our integrated reaction-diffusion system, the adjustable parameters were few,and the number of oscillations per trial and the numberof trials were large, thereby facilitating a detailed com-parison between theory and experiment. Having shownthe successes and limitations of the theorem using thegeneral methods of phase-reduction, it provides a frame-work for analyzing other networks. Most readily com-parable will be other networks with polygonal geometry.For larger-scale networks, it remains in principle possibleto numerically employ a similar methodology to predictspatiotemporal patterns by applying the H/K frameworkin complex networks, for example by using computer-assisted calculations [45]. However, for large numbers ofnodes and symmetries, it may become more practical andmeaningful to approximate the network by a continuumand use continuous symmetry groups (e.g., dense latticescould be approximated by planes, or polygons with alarge number of nodes by circles, with continuous rota-tions or translations as symmetries) [50–53]. Moreover,extensions of the theory will allow prediction of furtherinvariant manifolds, which do not arise directly from thesymmetry group of the network [10, 54, 55].These experiments raise for the first time the deep the-oretical question of how, in spite of this general consis-tency between symmetry based theories and experiment,that even small levels of heterogeneity have the potentialof crucially modifying the dynamics. Two examples arethat heterogeneity renders some symmetry-derived states no longer observable and the surprising phenomenonwhereby the sensitivity of states to heterogeneous per-turbations does not correlate with the strength and sizeof the basin of attraction of the state in the ideally sym-metric system. These results emphasize the importanceof assessing the robustness of symmetry-predicted resultsin the face of heterogeneity. This assessment is essentialto validate the application of the symmetry-based theoryto biological systems, as well as to guide the design ofchemical reaction-diffusion networks to be used in engi-neered applications, such as soft robotics.To date, we know of no theoretical framework address-ing the structural stability of H/K’s predictions to hetero-geneity that can explain our experimental and numericalobservations. Thus, this work encourages more theoreti-cal and experimental studies such as systematically intro-ducing symmetry breaking by controlling the degree towhich nodes and connections are distinct, so as to finelycharacterize the origin of heterogeneity-induced destabi-lization or vanishing of steady-states[28]. Our experimen-tal system, used here for the first time to test symmetrybased theories in reaction-diffusion networks, is ideallysuited for such studies[16]. Our results partially revealthe complex role of network structure on dynamics, butto articulate fully the engineering principles of networkdynamics it remains to elucidate how heterogeneity im-pacts performance. We hypothesize that similarly to thephase-locked and transient dynamics studied here, theimpact of network heterogeneity is partially symmetrygeneric and partially model specific.
ACKNOWLEDGEMENTS
We acknowledge financial support from NSF DMREF-1534890, the U. S. Army Research Laboratory and theU. S. Army Research Office under contract/ grant num-ber W911NF-16-1-0094, the microfluidics facility of theNSF MRSEC DMR-2011486, and the Swartz FoundationGrants 2017-6 and 2018-6.IH performed all experiments, data analysis, and simu-lations. MN guided all of the work. Experimental designby IH, MN, and SF. BC, CS, MN, JT, and IH contributedto understanding theoretical role of symmetry in system.MM helped fabricate microfluidic chips. IH, MN, JT,and SF wrote the manuscript.Simulations were performed using Brandeis Univer-sity’s High Performance Computing Cluster which is par-tially funded by DMR-MRSEC 2011486We acknowledge R´emi Boros, Youssef Fahmy, andAmanda Chisholm for their preliminary experiments on4 ring networks. We are grateful for Jan Engelbrecht andRennie Mirollo for their spirited discussions on dynamicalsystems theory.1
APPENDIX A: EXPERIMENTAL METHODSD. Network Fabrication
The microfluidic reaction-diffusion network was madeout of four adjacent reactors embedded in polydimethyl-siloxane (PDMS). The reactors are formed out of divotsin PDMS, forming effective buckets, which can be filledand then sealed all together by a piece of glass, form-ing a common lid. We manufactured these divots us-ing a soft lithographic process in which PDMS is curedwhile pressed against an inverse (positive) of the divotsmade out of a photoresist deposited onto a silicon wafer.This was performed as previously published [16], withthe exception of one adaptation described below. Thisgenerates a glass microscope slide coated with many re-actors organized into networks of four reactors, shown inFig. 1(b)(c) and Supplementary Material [29] Fig. 1 and2(a).The dimensions chosen for the network allow for ro-bust coupling of four nodes in ring topology. By ad-justing the sizes and distances between reactors wefound that rectangular reactor dimensions 62 µ m x 62 µ mx 30 µ m (L x W x H) with side-to-side distance 26 µ m resulted in strong coupling. The network reac-tors are organized in a 2 by 2 grid [Fig. 1(b)(c)] insuch a way that nearest neighbor reactors possess muchmore shared surface area relative next-nearest neigh-bors across the diagonal. This results in a ring-likeconnectivity, where coupling between nearest neighborsis stronger than across the diagonal. The rectangle ofBZ surrounding the network [Fig. 1 (c)] is forced into asteady-state, setting the concentration of chemicals sur-rounding the network. During each experiment we ob-serve nine or sixteen strongly coupled, individual net-works, separated from one another by controlled barriers[Supplementary Material [29] Fig. 1(b)(c)].The only alteration of the procedure in fabricat-ing the PDMS networks published [16] was to changethe way in which the PDMS was pressed and cured– instead of a 15kg lead brick applied for 12 hoursfollowed by baking in a 70C oven, we used a ther-mal press applying 90-113kg set at 70C for 2.5 hours.This was found to: a) reproducibly keep the size ofthe layer of PDMS underneath sample less than 2 µ m[Supplementary Material [29] Fig. 1(d)] and b) decreasethe probability that the silicon wafer breaks per use. E. Sample Holders
In a previous work the PDMS reactors had BZ sealedinside of them and were loaded into a microscope usingan acrylic plastic clamp [16]. This clamp did not controlthe temperature of the BZ. However, the frequency of BZoscillations depends on temperature [56].To maximize experimental reproducibility, wecreated a clamp that controlled sample tempera-
TABLE III. Final experimental chemical conditions in reac-tors: Chemical Molecular Formula Concentrationm m Sulfuric Acid H SO H O H FeN O S 3Tris(2,2’-bipyridyl)dichlororuthenium(II)hexahydrate C H Cl N Ru · O 1 . tures to within 0.1 ° C. The clamp’s temperatureis controlled through a thermistor that measuresthe temperature of the clamp nearby the sample[Supplementary Material [29] Fig. 2], 2 Peltier (TEC)devices [Supplementary Material [29] Fig. 2], and PIDfeedback between them mediated by an Arduino. Thesample is robustly driven to the clamp’s temperaturebecause the clamp possesses a large thermal mass rela-tive the sample and large thermal contact area with thesample [Supplementary Material [29] Fig. 2]. During alltrials samples were kept at 22 . ° C.We seal samples in the temperature-controlled clampexactly the same way as with the previous, plastic clamp[16], described in Appendix II C Protocol.
F. BZ Chemical Preparation
The BZ loaded into the microfluidic network is firstmixed outside the microfluidic device. A .24mL vol-ume of photo-sensitive BZ is prepared by sequentiallyadding equal 60 µ L volumes of Sulfuric acid, SodiumBromide, Malonic acid, Sodium Bromate, Ferroin thenTris(2,2’-bipyridyl)dichlororuthenium(II)hexahydrate toan Eppendorf tube, then mixing it with a Vortex mixer.Note that during the sequential pipetting of the chemi-cals, upon adding the Sodium Bromate, the solution con-verts from colorless to a vivid, transparent yellow for 15seconds before returning to a colorless state. The vol-umes output by the pipette used had a measured percentcoefficient of variance of 1 . G. Optics
We measured the chemical state of the reactorsthrough measuring their absorbance of green light. Fer-roin’s absorbance of green light changes drastically be-tween its oxidized and reduced state. The green light isfiltered to 515 ± .2Light perturbation, used to set boundary conditionsand initial conditions were set by projection of patternedblue light onto the sample, selectively exciting the pho-tocatalyst. As in previous works [16, 18, 23], the pat-terned blue light was periodically turned off and on ata high frequency to allow accurate measurement of ab-sorbance of Ferroin, with period 2 seconds and duty cycleof 50%. During some of the experiments the blue lightwas homogenized by directly replacing the sample with aCCD and using feedback between projected signals andthe measured values to minimize measured variability, asprevious published [19].Boundary conditions were applied by shining light onthe rectangle surrounding the network at an intensitythat completely inhibit oscillations in it.Initial conditions were set by applying light to the re-actors by inhibiting all reactors with light for 300-600seconds. Then, the light was turned off at different timesfrom each of the reactors, thus causing them to resumeoscillating at different times. The success rate of hittingtarget initial conditions far from Trot or Pronk was low.The light intensity of sample illumination was mea-sured by placing a power meter in the sample plane,the results are similar to in previous [23] work: In-tensity of blue light applied to boundaries: 0 . ± . − , Intensity of blue light applied to reactorsduring initial condition setting: 1 ± . − , In-tensity of blue applied light when projector blank/black:0 . ± . − , Intensity of 515nm green sampleillumination: ∼ − . Errors, in standard devia-tions, express variance in average illumination across thewhole sample field of view across all experiments, not thevariance across the field of view in individual experiments H. Protocol
The protocol for an experiment is as follows:1. PDMS chip, reentrant window, and O-ring[Supplementary Material [29] Fig. 2(a)] arecleaned with isopropyl alcohol, deionized water,and dried with compressed air. They are left underpetri dishes to prevent dust accumulation.2. A small batch of BZ solution is prepared as detailedearlier in Appendix II C. Solution is left in a darkchamber.3. The PDMS chip is plasma treated for 3 minutes at400mbar in ambient atmosphere.4. The BZ solution is then pipetted into the networksof interest in the PDMS chip as shown in depth insupplementary movie S7 of [16].5. Now, with the reentrant window placed ap-proximately above a feature of networks cov-ered by BZ, the reentrant window must besecured more firmly and precisely. While viewing the sample using a stereomicroscopewith green filtered transmission illumination, thethumbscrews [Supplementary Material [29] Fig. 2]are slowly turned, clamping the device. We al-ternated tightening them in a zig-zag pattern,with each tightening of a screw being roughlya 1/8 or less rotation. During this processany bubbles which are present in the reactorsshould decrease in size until they are invisible.Once all reactors are surrounded by dark out-lines Supplementary Material [29] Fig. 1(a), thereare no shearing distortions to the network, andthere are no bubbles, this process is halted.6. The clamp and the network with BZ sealed into itare then left in a dark, room temperature cham-ber until it has been 40 minute since the BZ wasinitially mixed in step 2, typically 20 minutes.7. The clamp is then loaded into theprojection illumination microscopeSupplementary Material [29] Fig. 3(c). Then,a MATLAB code with GUI is used toalign a projected pattern onto the sampleSupplementary Material [29] Fig. 3(b) and initiatetemperature control.8. Light is projected onto boundaries and sets ini-tial conditions of networks as described earlier inAppendix II C. Data is gathered for between 3000and 24000s, ∼
10 and 81 periods of oscillation ofeach reactor.9. In a few experiments a second attempt at settinginitial conditions was made.
APPENDIX B: PHASE-LOCKED CRITERIA
To identify phase-locked states in experiments we re-quire that ddt ( φ i − φ j ) is almost zero and is not acceler-ating, d dt ( φ i − φ j ) is also small. The algorithm used:1. Calculate the three phase differences versus time:( θ , θ , θ )2. Lowpass them to form: ( θ , θ , θ )3. Find the longest region in the ( θ , θ , θ ) timeseries when their velocities at below as threshold: | ddt θ | , | ddt θ | , | ddt θ | < . × − (cid:20) rads (cid:21) and the average acceleration is also below a thresh-old: ddt
13 ( | ddt θ | + | ddt θ | + | ddt θ | ) < × − (cid:20) rads (cid:21)
4. If the longest region is 5 or more periods of oscil-lation of all oscillators (1500 seconds), we considerthe experiment to be phase-locked.3
APPENDIX C: BEST FIT MODELI. Physical Models of BZ Microoscillators
The chemical concentration oscillations of an isolatedBZ microreactor are accurately modeled by the reactionkinetics derived for macroscopic reactors [19, 30]. De-noting concentrations, ¯ c = ( x, y, z, u ) where x =[HBrO ], y =[Br – ], z =[Oxidized catalyst], u =[Br ]: ddt ¯ c = ¯ R (¯ c ) (5)¯ R ( xyzu ) = R x ( x, y, z ) R y ( x, y, z, u ) R z ( x, z ) R u ( x, y, u ) R x ( x, y, z ) = k y − k xy − k x + k x ( c o − z )( c o − z + c min ) R y ( x, y, z ) = − k y + k u + k z − k xy − k x + k I ( c o − z )( b C b + 1) R z ( x, z ) = − ( k + k ) z + 2 k x ( c o − z )( c o − z + c min ) + k I ( c o − z )( b C b + 1) R u ( x, y, u ) = k y − k u + 2 k xy + k x The exchange of chemicals between adjacent BZ mi-croreactors is limited because the PDMS between them isapolar. Br is the only apolar intermediate of the BZ re-action which is certainly soluble in and diffusing throughthe PDMS between reactors. HBrO may also be solubleto a lesser degree [18, 24, 31]. Because of the short sepa-rations between reactors O (10 µm ), the diffusion of thesechemicals between reactors is assumed to be quasi-static[22, 28]. This results in a simple form of linear differencecoupling between a pair of BZ microreactors, where k isthe diffusive coupling rate of Br and k e is the ratio ofthe diffusive coupling rate of HBrO relative Br : ddt ¯ c = ¯ R (¯ c ) + µ (¯ c − ¯ c ) ddt ¯ c = ¯ R (¯ c ) + µ (¯ c − ¯ c ) µ = k e k k To simulate our results, we adapt the above, estab-lished model. Firstly, we model the ring-like topology byadding strong diffusive coupling between nearest neigh-bors and weaker coupling coupling across the diagonal.Secondly, we also permit a small degree of heterogene-ity in reaction rates within the reactors, representing:
TABLE IV. Simulation parameters known:Reagent concentrations:Description Value Unit a Bromate 288 m m m Malonicacid 400 m m c o Totalmetalioncatalyst 4 . m h Protons 160 m m b Bromomalonicacid 0 . ∗ m m m Reagent rates and relevant constants:Value Unit k × h m − s − k h a s − k × m − s − k ha s − k × h m − s − k
10 s − k m s − k . m s − k b s − k . m s − k r × m − s − k red × m − s − k I − b C . m c min p k r ( k + k ) c o /k red m i) differences in reagent concentrations, ii) differences inpowers of light hitting photo-sensitive reactors, and iii)slightly different boundary conditions. The model is thena function of: a) k , the diffusive coupling rate of Br , b) k e , the ratio of the coupling rate of HBrO relative Br c) f , the ratio of diagonal to nearest neighbor coupling,d) differences in local reactor reaction rates: ddt ¯ c i = ¯ R o (¯ c i ) + (cid:15) ¯ r i (¯ c i ) + X j = i A ij µ (¯ c j − ¯ c i ) (6) A ij = f
11 0 1 ff f ; µ = k e k k ;To ease comparison to experiments and qualitativelyunderstand the role of free parameters, we use establishedmethods [36, 38, 39], described later in this section, to re-duce the model above Eqn. 6 to the phase model Eqn. 4.The reduced model Eqn. 4 predicts the phase dynamicsof the reactors relative one another as a function of theparameters: a) k which scales the interaction function H , b) k e which determines the shape of the H functionFig. 6(a), c) f which changes the adjacency matrix, d)4differences in local reactors which change intrinsic fre-quencies ω i .We then construct a 3D model of the evolution ofphase-differences between adjacent nodes in the network.Letting θ ij ≡ φ i − φ j and ∆ ω ij ≡ ω i − ω j ddt θ ij = ∆ ω ij + k (cid:20) X k = i A ik H ( φ k − φ i ) − X k = j A jk H ( φ k − φ j ) (cid:21) (7)(8) ddt ¯ θ = ddt θ ddt θ ddt θ = ¯Ψ(¯ θ ) (9)This gives us a model which is expected to qualitativelycapture the dynamics of the system and allow us to betterunderstand it. Since the parameters of the model, f , k and so on, cannot be independently measured, the modelmust be empirically fit to experimental data. However,the best fit model reveals important, nontrivial insightsinto our particular symmetric 4 node network and BZmicrooscillator networks in general.The H function is shown in Fig. 6(a). We ran simula-tions with the H function saved as a Chebychev function[57] using MATLAB’s ODE45 with relative and absolutetolerances of 1 × − . J. Fitting to Experiments
The model Eqn. 9 was fitted to each experimentaltime series of phase differences versus time [Fig. 6(b)]using nonlinear regression. The nonlinear regressionwas performed on each individual experiment by con-straining a simulation to start from an experimentalinitial condition, then optimizing the simulation’s pa-rameters to reduce the squared error between its tra-jectory and the experiment’s trajectory using matlab2019b’s surrogate optimization. Specifically, the ini-tial condition of phase difference of an experiment wasset to be the phase differences when the reactors meetthe similar frequency threshold defined in AppendixII C, shown in Supplementary Material [29] Fig. 4(a)(b).The final point in an experimental trajectory usedin a fitting was half way between when the phase-locked condition was met and when it was lost[Supplementary Material [29] Fig. 4(a)(c)] or the end ofthe experiment if it did not unlock.The identified best fit parameters were: ∆ ω ij = 0 , k =1 . × − ± . × − [s − ], f = 0[1] and k e = . ω ij obeyeda Laplacian distribution ρ (∆ ω ij ) = b exp( − | ∆ ω ij − µ | b ) TABLE V. Simulation parameters fitted: In ’Fit values’ and’Values used in theory’ a single number represents the num-ber fit or used in simulations. If in a fit or simulation valueswere randomly distributed, the form of the distribution is de-scribed by ’G’ or ’L’. ’G( µ, σ )’ represents a Gaussian proba-bility density function with mean µ and standard deviation σ .A ’L( µ, b )’ represents a Laplacian probability density functionwith mean µ and rate parameter b .Parameter Fit values Values used intheory UnitCoupling: k G(1 . × − , . × − ) 2 × − s − k e f ω ij L(0 , π × − ) L(0 , π × − ) rad s − [Fig. 6(c)]. The Laplacian distribution of intrinsic fre-quencies has a mean, µ , of 0 and rate parameter, b , of2 π × − [rad s − ]. Further, the fitting required someexcitatory coupling k e = . f = 0[1]. The best fit values and the values used insimulations shown in Table V.For a discussion of the values of these parameters mea-sured in previous works and their possible physical rele-vance, please see Supplementary Material [29] Sec. II. K. Best Fit Simulations with Heterogeneity
To determine the impact of experimentally realisticheterogeneity on the model, we ran simulations of Eqn. 4with the experimental best fit parameters in Table V.Specifically, heterogeneity in intrinsic frequencies in sim-ulations ∆ ω ij were drawn from the distribution which fitsexperiments, a Laplacian distribution defined in TableV, with mean 0 and rate parameter 2 π × − [rad s − ],while all other parameters are the exact, constant valueenumerated in the third row of Table V.In running simulations from a dense set of initial con-ditions, each initial condition had 7 simulations initial-ized from it with independent resamplings of frequencyheterogeneity. We thus could observe the impact of het-erogeneity throughout state-space by sampling the distri-bution of heterogeneity in all regions of state-space. Theresult of 34,713 such simulations are shown in Fig. 5(b). L. Computing Phase Model Reduction
We use established methods [36, 38, 39] to reduce thereaction-diffusion model of our 4 reactor network Eqn. 6to the phase model Eqn. 4. In this framework, we firstdetermine the phase-dependent phase shift of an uncou-pled BZ reactor induced by sudden, small additions of itschemical species. We then compute the fluxes of chemicalspecies between a pair of diffusively coupled reactors as5
ExperimentSimulation θ θ θ HK steady state (b) Fit to an experiment (c) Laplace distribution fit intrinsic frequencies H Br H HBrO /10H
Best fit2 2 H ( χ ) [ r ad ] χ[rad] (a) Coupled BZ reactor interaction functions Heterogeneity ofbest fit to experiment C oun t Time[ks] [ r ad ] P ha s e d i ff e r en c e Exponential pdf fit
FIG. 6. (a)
The interaction functions of coupled BZ oscillators. H Br is the interaction functions of purely Br coupledreactors. H HBrO is the interaction functions of purely HBrO coupled reactor. In the model of our 4 node network Eqn. 9,when k e is not zero the interaction function is H ( χ ) = H Br ( χ )+ k e H HBrO ( χ ). Since H HBrO is largely 0 except for χ between 0and π/ k e adds a bump of phase advance near 0, without affecting the otherwise phase delaying dynamics due to Br . H Best fit corresponds to the interaction function which best fits experiments k e = 0 . (b) Example experimental trajectory convergingto steady-state (circles), and best fit phase model (solid line). The fit model required slight heterogeneity ∆ ω ij = 0. Thesteady-state phase difference ( π, π, π ) of the phase model, without best fit heterogeneity, at dashed horizontal line. This is theH/K predicted point invariant manifold nearest to the experimental trajectory. (c) Distribution of best fit nondimensionalizedintrinsic frequency differences. They are nondimensionalized by dividing them by best fit coupling rate k times the amplitudeof the interaction function H , H max . a function of their relative phases. By combining these,we determine the rate of phase shift of a reactor as afunction of its phase difference with a neighbor.If chemicals are added to a BZ reactor on its limit cycle,the time until it next spikes will be changed. These timeshifts, or phase shifts, depend on amount of chemicalsadded as well as the time along a period of oscillation, orphase, they are added at. We quantify the phase shifts ofan uncoupled reactor to a perturbation at a given phasein terms of phase response curves (PRC). The PRC, Q ,is a set of curves of phase shift, ∆ φ , as a function of achemical perturbation’s chemical concentration, ∆¯ c , andthe phase, φ , at which it is added at. Specifically, itallows calculation of ∆ φ = Q ( φ ) . ∆¯ c . the case of the 4Dmodel of a BZ reaction 5, Q : R S is a set of 4curves with units of phase per chemical concentrations[rad m − ]. We compute Q in the limit of infinitesimalperturbations | ∆¯ c | → Q as the total derivative Q ( φ ) = ∇ ¯ c φ .Noting that rate of change of phase of a reactor is ddt φ i and is a sum of constant term, its intrinsic frequency ω i , and a time-varying term due to its coupling with itsneighbors depending on their relative phases: ddt φ i = ω i + X j A ij F ( φ i , φ j )Using the PRC we can then compute the rate of change of phase of a reactor F [rad s − ] with respect to an incident,relative phase dependent flux due to coupling ¯ g ( φ i , φ j )[ m s − ] using the PRC and the chain rule: F ( φ i , φ j ) = Q ( φ )¯ g ( φ i , φ j ) ≈ ∇ ¯ c φ ¯ g ( φ i , φ j ) (10)Letting the unperturbed limit cycle of a BZ reactor beparametrized ¯ c LC ( φ ) and denoting the neighboring re-actor j , we see from our diffusive coupling Eqn. 6 thatthe specific form of phase dependent flux neighbors ex-perience ¯ g ( φ i , φ j ) = µ [¯ c LC ( φ j ) − ¯ c LC ( φ i )] and it followsthat: F ( φ i , φ j ) = Q ( φ ) µ [¯ c LC ( φ j ) − ¯ c LC ( φ i )] (11)Note that F is a linear combination of flux due to Br , u , and HBrO , x : F ( φ i , φ j ) = k [ Q u ( u LC ( φ j ) − u LC ( φ i ))+ k e Q x ( x LC ( φ j ) − x LC ( φ i ))] (12)Since the change of relative phase differences is smallduring a period of reactor oscillation, we compute F ( φ j , φ i ) averaged over a cycle [36, 38, 39]. Doing sotransforms it into a function of relative phase difference H ( φ j − φ i ). Since F is scaled by k , we define H such6that it must be used by explicit scaling by k , the diffu-sive coupling rate of Br [s − ] in eqn 4: H ( φ j − φ i ) ≡ k − (2 π ) − Z π F ( α, α + φ j − φ i ) dα (13)We express the interaction function H as the sum of twodistinct, separately calculable terms, because F is a linearcombination of two functions, one for Br and one forHBrO : H ( φ j − φ i ) = H Br ( φ j − φ i ) + k e H HBrO ( φ j − φ i ) (14) APPENDIX D: METRIC OF DISTANCE
To measure distances between two points in the state-space of the 3D phase difference dynamics, both experi-mentally and in simulations, we found a surprising func-tion d (¯ θ, ¯ θ ) was required. For a given pair of points ¯ θ and ¯ θ d is calculated by the following algorithm:1. Consider two points in the state-space: ¯ θ =( θ , θ , θ ) and ¯ θ = ( θ , θ , θ )2. Compute phase difference of fourth edge,which is completely determined by the otherthree ¯ θ f = ( θ , θ , θ , θ + θ + θ ) =( θ , θ , θ , θ ) and ¯ θ f = ( θ , θ , θ , θ )3. Define a vector of phase difference between stateswith ∠ being complex, or phasor, angle θ diff j = ∠ exp( i ∗ ( θ f j − θ f j ))4. Let the distance between ¯ θ and ¯ θ be the Euclideannorm of the 4D phase difference vector d (¯ θ, ¯ θ ) = | ¯ θ diff | APPENDIX E: COMPUTING TRANSVERSELYAPUNOV EXPONENTS OF INVARIANTMANIFOLDS
The maximum transverse Lyapunov exponent (MTLE)of the higher order invariant manifolds describe whethertrajectories collapse or diverge from them. If a higher or-der invariant manifold has a negative MTLE, trajectorieswill converge towards it [43]. In our particular system,we found the MTLE were often easy to calculate and easyto qualitatively understand.The MTLE of the invariant manifolds can be computedat the phase model level using the well known methodof linearizing the dynamics Eqn. 9 about the invariantmanifolds [43–45, 58, 59]. We explicitly list how we per-formed this procedure in our case in which the invariantmanifolds are linear hyperplanes.We begin by considering whether a point ¯ θ perturbedoff an invariant manifold, displaced by (cid:15) ¯ δ from point ¯ θ ,converges back to or diverges from the manifold [Fig. 7]. Dynamics near a linear invariant manifold:
Trajectory on invariant manifoldDisplacementLinear invariant manifoldTrajectory nearinvariant manifoldInvariant manifoldtangentInvariant manifoldnormal
FIG. 7. Diagram of flows nearby a linear invariant man-ifold. Two example trajectories ¯ θ and ¯ θ , which are flowsin the state-space of Eqn. 9. Note that because ¯ θ is on aninvariant manifold, it must continue flowing tangent the in-variant manifold indefinitely. However, the trajectory nearthe invariant manifold ¯ θ is not limited in this way. We canconsider the displacement between these two trajectories as afunction of time (cid:15) ¯ δ ( t ). In the case shown, the displacementis shrinking, indicating that the nonlinear flow about the lin-ear invariant manifold near ¯ θ is converging towards it. Toquantify this behavior more broadly, we utilize he fact thenormal component of the displacement will take on a form¯ n T (cid:15) ¯ δ ( t ) ∝ exp( λt ). The exponential scale λ , called the trans-verse Lyapunov exponent, when negative indicates attraction,when positive repulsion. The dynamics of the perturbed point take on a simpleform: ddt ¯ θ = ¯Ψ(¯ θ ) = ¯Ψ(¯ θ + (cid:15) ¯ δ )= ¯Ψ(¯ θ ) + (cid:15) ∇ ¯ x ¯Ψ(¯ x ) | ¯ x =¯ θ ¯ δ + O ( (cid:15) ) (15)It is convenient to consider the displacement between theperturbed solution and the solution on the invariant man-ifold: ddt ¯ δ = ddt ¯ θ − ddt ¯ θ = ¯Ψ(¯ θ ) − ¯Ψ(¯ θ )= (cid:15) ∇ ¯ x ¯Ψ(¯ x ) | ¯ x =¯ θ ¯ δ + O ( (cid:15) ) (16)We can determine if the component of the perturbationalong a given direction normal to the invariant manifold,¯ n , grows or decays of the displacement off of the invari-ant manifold. Letting J (¯ x ) ≡ ∇ ¯ x ¯Ψ(¯ x ) | ¯ x =¯ x , the rate ofchange of the normal distance, ¯ n T ¯ δ is given by: ddt (¯ n T ¯ δ ) ≈ ¯ n T (cid:15)J (¯ θ )¯ δ (17)7To precisely determine the exponential timescale withwhich the normal components of the perturbed trajectorygrows or shrinks, we use a unitary transformation ma-trix P to projects a vector of phase differences from thecanonical basis into a basis aligned with a given invari-ant manifold: Letting ¯ t i be tangent the invariant mani-fold and ¯ n i be perpendicular, P = [¯ t ... ¯ n ... ] T s.t. if ¯ θ = c t ¯ t + c t ¯ t ... + c n ¯ n + .. . then P. ¯ θ = ( c t , c t , ...c n ... ).Crucially, points on the invariant manifold have no nor-mal components ∀ j c nj = 0 thus P. ¯ θ = ( c t , c t , ... , .. ).Similarly, on the invariant manifold the velocity field P. ¯Ψ(¯ θ ) = ddt ( c t , ...c n , ... ) must have all its normal com-ponent be 0, ∀ j ddt c nj = 0, as by definition trajectoriescannot flow out of, or transverse, an invariant manifoldFig. 7.We can now consider the growth or decay of perturba-tions in or out of a manifold. Let ¯ ξ = P. ¯ δ : ddt ¯ ξ = ddt ( P. ¯ δ ) = P J ¯ δ = P JP − ¯ ξ (18)Letting J = P JP − : ddt ¯ ξ = J ¯ ξ - a linear time-varyingequation. Since the dynamics of the system ¯ ξ normal tothe invariant manifold are all 0, independent on the loca-tion along the manifold spanned by the tangent, it followsthat the block dynamics in J corresponding to normalcomponents will be independent of those for tangent com-ponents [43]. The maximum real eigenvalue of this blockof normal dynamics gives an approximate bound on theexponential growth or decay of the trajectory’s distancenormal an invariant manifold and is called the maximumtransverse Lyapunov exponent. M. Algorithm
The algorithm, detailed at length above, is as follows:1. Choose an invariant manifold M , of dimension k ,in an n node network.2. Determine k orthonormal vectors which span M and label them the tangent vectors T . Determine aset N of n − − k vectors orthonormal one anotherand T . We use the Gramm-Schmitt procedure.3. Compute a unitary transformation matrix P , withthe first columns composed of invariant manifoldtangents, then followed by normals.4. Compute the Jacobian of your nonlinear flow atpoints on the invariant manifold J (¯ θ ), written ex-plicitly as a function of location on the invariantmanifold ¯ θ .5. Transform the Jacobian into its tangent and normalcomponents using P via P.J (¯ θ ) .P − = J (¯ θ ).6. Extract the block of J that contains the decou-pled transverse dynamics - the columns and rowscorresponding to normal components.7. Compute the maximum real eigenvalue of the nor-mal block λ (¯ θ ), which is the max TransverseLyapunov Exponent An example of executing this algorithm for the( D p , D p ) invariant manifold of the 4 ring is inSupplementary Material [29] Sec. IIIA. Generic expres-sions of block of J for all invariant manifolds of a broadclass of 4 ring networks are computed in terms of firstderivatives of H in Supplementary Material [29] TableI. A comparison of the MTLE of a ring of 4 Ku-ramoto oscillators to our system is presented inSupplementary Material [29] Fig. 8. [1] S. H. Strogatz, “Exploring complex networks,” Nature , 268–276 (2001).[2] Adilson E. Motter, Seth A. Myers, Marian Anghel, andTakashi Nishikawa, “Spontaneous synchrony in power-grid networks,” Nature Physics , 191–197 (2013).[3] Martin Golubitsky, Ian Stewart, Pletro Luclano Buono,and J. J. Collins, “Symmetry in locomotor central pat-tern generators and animal gaits,” Nature , 693–695(1999).[4] Atsuko Takamatsu, Reiko Tanaka, Hiroyasu Yamada,Toshiyuki Nakagaki, Teruo Fujii, and Isao Endo, “Spa-tiotemporal symmetry in rings of coupled biological os-cillators of physarum plasmodial slime mold,” PhysicalReview Letters , 078102 (2001).[5] Visarath In, Andy Kho, Joseph D. Neff, Antonio Pala-cios, Patrick Longhini, and Brian K. Meadows, “Exper-imental observation of multifrequency patterns in arraysof coupled nonlinear oscillators,” Physical Review Letters , 244101 (2003).[6] Matthew H. Matheny, Jeffrey Emenheiser, Warren Fon,Airlie Chapman, Anastasiya Salova, Martin Rohden,Jarvis Li, Mathias Hudoba de Badyn, M´arton P´osfai,Leonardo Duenas-Osorio, Mehran Mesbahi, James P.Crutchfield, M. C. Cross, Raissa M. D’Souza, andMichael L. Roukes, “Exotic states in a simple net-work of nanoelectromechanical oscillators,” Science ,eaav7932 (2019).[7] Ian Stewart and Martin Golubitsky, “Symmetry methodsin mathematical biology,” S˜ao Paulo J. Math. Sci. , 1–36 (2015).[8] P. Ashwin and J. W. Swift, “The dynamics of n weaklycoupled identical oscillators,” Journal of Nonlinear Sci-ence , 69–108 (1992).[9] Martin Golubitsky and Ian Stewart, The Symmetry Per-spective From Equilibrium to Chaos in Phase Space andPhysical Space , 1st ed. (Birkhauser Verlag, Basel · Boston · Berlin, 2000) pp. 59–79.[10] Martin Golubitsky and Ian Stewart, “Nonlinear dynam-ics of networks: The groupoid formalism,” Bulletin of theAmerican Mathematical Society , 305–364 (2006).[11] Martin Golubitsky and Ian Stewart, “Rigid patterns ofsynchrony for equilibria and periodic cycles in networkdynamics,” Chaos , 094803 (2016).[12] Einat Couzin-Fuchs, Tim Kiemel, Omer Gal, Amir Ayali,and Philip Holmes, “Intersegmental coupling and recov-ery from perturbations in freely running cockroaches,”Journal of Experimental Biology , 285—-297 (2015).[13] Calvin Zhang, Robert D Guy, Brian Mulloney, Qing-hai Zhang, and Timothy J Lewis, “Neural mechanismof optimal limb coordination in crustacean swimming,”Proceedings of the National Academy of Sciences ,13840–13845 (2014).[14] N. Kopell and G. B. Ermentrout, “Coupled oscillatorsand the design of central pattern generators,” Mathe-matical Biosciences , 87–109 (1988).[15] Martin Golubitsky and Ian Stewart, The Symmetry Per-spective From Equilibrium to Chaos in Phase Space andPhysical Space , 1st ed. (Birkhauser Verlag, Basel · Boston · Berlin, 2000) pp. 1–337.[16] Thomas Litschel, Michael M. Norton, Vardges Tserun-yan, and Seth Fraden, “Engineering reaction-diffusionnetworks with properties of neural tissue,” Lab on a Chip , 714–722 (2018), arXiv:arXiv:1711.00444.[17] Alan Turing, “The chemical basis of morphogenesis,”Philosophical Transactions of the Royal Society of Lon-don. Series B, Biological Sciences , 37–72 (1952).[18] Nathan Tompkins, Ning Li, Camille Girabawe, MichaelHeymann, G Bard Ermentrout, Irving R Epstein, andSeth Fraden, “Testing Turing’s theory of morphogenesisin chemical cells,” Proceedings of the National Academyof Sciences , 4397–4402 (2014).[19] James Sheehy, Ian Hunter, Maria Eleni Moustaka, S. AliAghvami, Youssef Fahmy, and Seth Fraden, “Impactof pdms-based microfluidics on belousov–zhabotinskychemical oscillators,” The Journal of Physical ChemistryB (2020), 10.1021/acs.jpcb.0c08422.[20] Atsuko Takamatsu, “Spontaneous switching among mul-tiple spatio-temporal patterns in three-oscillator systemsconstructed with oscillatory cells of true slime mold,”Physica D: Nonlinear Phenomena , 180–188 (2006).[21] Alexandra M. Tayar, Eyal Karzbrun, Vincent Noireaux,and Roy H. Bar-Ziv, “Synchrony and pattern formationof coupled genetic oscillators on a chip of artificial cells,”Proceedings of the National Academy of Sciences of theUnited States of America , 11609–11614 (2017).[22] Masahiro Toiya, Hector O. Gonz´alez-Ochoa, Vladimir K.Vanag, Seth Fraden, and Irving R. Epstein, “Synchro-nization of chemical micro-oscillators,” Journal of Phys-ical Chemistry Letters , 1241–1246 (2010).[23] Jorge Delgado, Ning Li, Marcin Leda, Hector O.Gonz´alez-Ochoa, Seth Fraden, and Irving R. Epstein,“Coupled oscillations in a 1D emulsion of Belousov-Zhabotinsky droplets,” Soft Matter , 3155–3167 (2011).[24] Ning Li, Jorge Delgado, Hector O Gonz´alez-Ochoa,Irving R Epstein, and Seth Fraden, “Combined ex-citatory and inhibitory coupling in a 1-D array ofBelousov–Zhabotinsky droplets,” Physical ChemistryChemical Physics , 10965–10978 (2014).[25] Nathan Tompkins, Matthew Carl Cambria, Adam L.Wang, Michael Heymann, and Seth Fraden, “Cre- ation and perturbation of planar networks ofchemical oscillators,” Chaos: An InterdisciplinaryJournal of Nonlinear Science , 064611 (2015),https://doi.org/10.1063/1.4922056.[26] Ning Li, Nathan Tompkins, Hector Gonzalez-Ochoa,and Seth Fraden, “Tunable diffusive lateral inhibition inchemical cells,” European Physical Journal E , 1–12(2015), arXiv:15334406.[27] A. L. Wang, J. M. Gold, N. Tompkins, M. Heymann,K. I. Harrington, and S. Fraden, “Configurable NORgate arrays from Belousov-Zhabotinsky micro-droplets,”European Physical Journal: Special Topics , 211–227(2016).[28] Michael M. Norton, Nathan Tompkins, Baptiste Blanc,Matthew Carl Cambria, Jesse Held, and Seth Fraden,“Dynamics of Reaction-Diffusion Oscillators in Star andother Networks with Cyclic Symmetries Exhibiting Mul-tiple Clusters,” Physical Review Letters , 148301(2019).[29] “See supplemental material at [url will be inserted bypublisher] for movies and details concerning experimentalmethods, numerical simulations, and theory.”.[30] Vladimir K. Vanag and Irving R. Epstein, “A model forjumping and bubble waves in the Belousov-Zhabotinsky-aerosol OT system,” Journal of Chemical Physics ,104512 (2009).[31] Kristian Torbensen, Federico Rossi, Sandra Ristori, andAli Abou-Hassan, “Chemical communication and dy-namics of droplet emulsions in networks of Belousov-Zhabotinsky micro-oscillators produced by microflu-idics,” Lab on a Chip , 1179–1189 (2017).[32] Ivan S. Proskurkin and Vladimir K. Vanag, “Dynamicsof a 1D array of inhibitory coupled chemical oscillatorsin microdroplets with global negative feedback,” PhysicalChemistry Chemical Physics , 16126–16137 (2018).[33] Kristian Torbensen, Sandra Ristori, Federico Rossi, andAli Abou-Hassan, “Tuning the Chemical Communica-tion of Oscillating Microdroplets by Means of MembraneComposition,” Journal of Physical Chemistry C ,13256–13264 (2017).[34] Vladimir K. Vanag and Irving R. Epstein, “Excitatoryand inhibitory coupling in a one-dimensional array ofBelousov-Zhabotinsky micro-oscillators: Theory,” Phys-ical Review E - Statistical, Nonlinear, and Soft MatterPhysics , 066209 (2011).[35] Arthur T. Winfree, “Biological rhythms and the behaviorof populations of coupled oscillators,” Journal of Theo-retical Biology , 15–42 (1967).[36] George Bard Ermentrout and David H Terman, Inter-disciplinary Applied Mathematics , 1st ed., edited by S.S.Antman, J.E. Marsden, L. Sirovich, and S. Wiggins,Vol. 8 (Springer, New York Dordrecht Heidelberg Lon-don, 2009) pp. 171–216.[37] Yoshiki Kuramoto,
Chemical Oscillations, Waves, andTurbulence , 1st ed. (Springer-Verlag, Berlin, 1984) p. 156.[38] Michael A. Schwemmer and Timothy J. Lewis, “The the-ory of weakly coupled oscillators,” in
Phase ResponseCurves in Neuroscience: Theory, Experiment, and Anal-ysis (2012) pp. 3–30.[39] Dan Wilson and Bard Ermentrout, “Augmented phasereduction of (Not So) weakly perturbed coupled oscilla-tors,” SIAM Review , 277–315 (2019).[40] Bharat Monga, Dan Wilson, Tim Matchen, and JeffMoehlis, “Phase reduction and phase-based optimal con- trol for biological systems: a tutorial,” Biological Cyber-netics , 11–46 (2019).[41] Steven Strogatz, Computers in Physics , 1st ed. (PerseusBooks Publishing L.L.C., New York, 1994) p. 154.[42] Krishna Pusuluri, Sunitha Basodi, and AndreyShilnikov, “Computational exposition of multistablerhythms in 4-cell neural circuits,” Communications inNonlinear Science and Numerical Simulation , 105139(2020).[43] Peter Ashwin, Jorge Buescu, and Ian Stewart, “Fromattractor to chaotic saddle: A tale of transverse instabil-ity,” Nonlinearity , 703–737 (1996).[44] Louis M. Pecora and Thomas L. Carroll, “Master stabil-ity functions for synchronized coupled systems,” PhysicalReview Letters , 2109–2112 (1998).[45] Louis M. Pecora, Francesco Sorrentino, Aaron M. Hager-strom, Thomas E. Murphy, and Rajarshi Roy, “Clustersynchronization and isolated desynchronization in com-plex networks with symmetries,” Nature Communica-tions , 4079 (2014).[46] Atsuko Takamatsu, Teruo Fujii, and Isao Endo, “Timedelay effect in a living coupled oscillator system with theplasmodium of physarum polycephalum,” Physical Re-view Letters , 2026–2029 (2000).[47] Per Sebastian Skardal, Dane Taylor, and Jie Sun, “Opti-mal synchronization of complex networks,” Physical Re-view Letters , 144101 (2014), arXiv:1402.7337.[48] Igor Belykh and Martin Hasler, “Mesoscale and clustersof synchrony in networks of bursting neurons,” Chaos ,016106 (2011).[49] Abu Bakar Siddique, Louis Pecora, Joseph D. Hart,and Francesco Sorrentino, “Symmetry- and input-clustersynchronization in networks,” Phys. Rev. E , 042217(2018).[50] I.B. Vivancos, P. Chossat, and I. Melbourne, “New plan-forms in systems of partial differential equations with Eu-clidean symmetry,” Archive for Rational Mechanics andAnalysis , 199–224 (1995). [51] William H. Bosking, Ying Zhang, Brett Schofield, andDavid Fitzpatrick, “Orientation selectivity and the ar-rangement of horizontal connections in tree shrew striatecortex,” Journal of Neuroscience , 2112–2127 (1997).[52] Paul C Bressloff, Jack D Cowan, Martin Golubitsky, andPeter J Thomas, “Scalar and pseudoscalar bifurcationsmotivated by pattern formation on the visual cortex,”Nonlinearity , 739–775 (2001).[53] P. C. Bressloff, J. D. Cowan, M. Golubitsky, P. J.Thomas, and M. C. Wiener, “Geometric visual hallu-cinations, Euclidean symmetry and the functional archi-tecture of striate cortex,” Philosophical Transactions ofthe Royal Society B: Biological Sciences , 299–330(2001).[54] Michael T. Schaub, Neave O’Clery, Yazan N. Billeh,Jean-Charles Delvenne, Renaud Lambiotte, and Mauri-cio Barahona, “Graph partitions and cluster synchroniza-tion in networks of oscillators,” Chaos , 94821 (2016),arXiv:1608.04283.[55] Anastasiya Salova and Raissa M. D’Souza, “Decoupledsynchronized states in networks of linearly coupled limitcycle oscillators,” Physical Review Research , 43261(2020), arXiv:2006.06163.[56] Tam´as B´ans´agi, Marcin Leda, Masahiro Toiya, Anatol M.Zhabotinsky, and Irving R. Epstein, “High-frequency os-cillations in the Belousov - Zhabotinsky reaction,” Jour-nal of Physical Chemistry A , 5644–5648 (2009).[57] Tobin A Driscoll, Nicholas Hale, and Lloyd N Trefethen, Chebfun Guide , 5th ed. (Pafnuty Publications, Oxford,2014).[58] V. N. Belykh, I. V. Belykh, and K. V. Nelvidin, “Spa-tiotemporal synchronization in lattices of locally coupledchaotic oscillators,” in
Mathematics and Computers inSimulation (2002).[59] Giovanni Russo and Jean Jacques E. Slotine, “Symme-tries, stability, and control in nonlinear systems and net-works,” Physical Review E - Statistical, Nonlinear, andSoft Matter Physics , 041929 (2011), arXiv:1011.0188. he Symmetry Basis of Pattern Formation in Reaction-Diffusion Networks Electronic Supplement
Ian Hunter, ∗ Michael M. Norton, ∗ Bolun Chen,
3, 4
Chris Simonetti, Maria Eleni Moustaka, Jonathan Touboul,
3, 5 and Seth Fraden † Brandeis University Physics, Waltham MA, 02453 USA Center for Neural Engineering, Department of Engineering Science and Mechanics,The Pennsylvania State University, University Park, PA 16802 USA Brandeis University Volen National Center for Complex Systems, Waltham MA, 02453 USA Department of Physics Boston University, Boston MA, 02215 USA Brandeis University Mathematics Department (Dated: March 2, 2021) ∗ These two authors contributed equally † To whom correspondence should be addressed:[email protected] a r X i v : . [ n li n . PS ] F e b I. MOVIES • Movie S1: Video of an experimental network converging to the Trot attractor: Initially three reactors flash at asimilar times, with the fourth flashing out of phase relative to them. As time goes on the system converges to anearly ideal realization of the Trot state. The raw video of the experiment is shown alongside a space-time plotof the video and the calculated phase differences between reactors 1 and 2, 2 and 3, and 3 and 4. The phasedifferences are denoted θ ij ≡ φ i − φ j . Link to video: https://tinyurl.com/3te6ex37 • Movie S2: Video of an experimental network in the Pace attractor: The experiment starts in a Pace pattern andproceeds to stay in it during the rest of the experiment. The raw video of the experiment is shown alongside aspace-time plot of the video and the calculated phase differences between reactors 1 and 2, 2 and 3, and 3 and4. Link to video: https://tinyurl.com/2d2p4d9u • Movie S3: Video of an experimental network in the Pronk attractor: The experiment starts in a Pronk stateand proceeds to stay in it the during rest of the experiment. The raw video of the experiment is shown alongsidea space-time plot of the video and the calculated phase differences between reactors 1 and 2, 2 and 3, and 3 and4. Link to video: https://tinyurl.com/apku3xcx • Movie S4: Video of an experimental network in the Other attractor: The experiment starts and stays in a statein which flashing is initiated at the top left, and proceeds to the bottom right reactor in a Z-pattern. Sincethis does not resemble any H/K state, it is categorized as Other. The raw video of the experiment is shownalongside a space-time plot of the video and the calculated phase differences between reactors 1 and 2, 2 and 3,and 3 and 4. Link to video: https://tinyurl.com/dupdjf9k • Movie S5: Video of rotation about 3D plot in manuscript Fig. 3(a) which shows the state-space of the theoreticalmodel. The video starts from the perspective of the Fig. 3A, then rotates about it to help readers understandthe 3D structure of the state-space. Link to video: https://tinyurl.com/63zsevcu • Movie S6: Video of rotation about a 3D plot of the same state-space as Fig. 3(a) except with experimental data.Two panels are shown: the left being the trajectories that converged Trot, Bound, and Pronk and the rightbeing the trajectories that converged Other. The left panel, has clear similarity to the simulated state-space inMovie S5. Link to video: https://tinyurl.com/c7n7tkmj • Movie S7: Video of rotation about 3D plot in manuscript Fig. 3(b) which shows basins of attraction of thetheoretical model and experiments. The video starts from the perspective of the Fig. 3(b), then rotates aboutit to help readers understand the 3D structure of the basins of attraction. Link to video: https://tinyurl.com/ye6thfme • Movie S8: Annotated video of experiment converging from Pronk to Pace Fig. 4(a). The raw video of theexperiment is shown alongside a space-time plot of the video and evolution of the phase differences betweenreactors moving through the state-space of the network. As during this experiment the 2nd and 3rd reactoroscillate in-phase, only the 2D slice of the state-space in which the trajectory is confined to, θ = 0, is shown.Within this view of state-space the 1D H/K predicted invariant manifold ( D s , D s ) is shown and is colored byits maximum transverse lyapunov exponent, as in the manuscript Fig. 4(a)(c). The experimental trajectorynoisily moves about Pronk before transitioning to Pace directly along the invariant manifold. Link to video: https://tinyurl.com/tsyw6xk3 • Movie S9: Annotated video of experiment converging from Pronk to Trot Fig. 4(b). The raw video of theexperiment is shown alongside a space-time plot of the video and evolution of the phase differences betweenreactors moving through the state-space of the network. As during this experiment the 1st and 3rd reactoroscillate in-phase, only the 2D slice of the state-space in which the trajectory is confined to, θ = − θ where φ = φ , is shown. Within this view of state-space the 2D H/K predicted invariant manifold ( D p , D p ) isshown and is colored by its maximum transverse lyapunov exponent, as in the manuscript Fig. 4(b)(c). Theexperimental trajectory starts near Pronk before efficiently transitioning to Trot, directly following the in-planevelocity field of the 2D the invariant manifold. Link to video: https://tinyurl.com/354k63ma • Movie S10: Video of rotation about 3D plot in manuscript Fig. 4(c) which shows theoretical stability of in-variant manifolds. The video starts from the perspective of the Fig. 4(c), then rotates about it to help readersunderstand the structure of the invariant manifold stability. Link to video: https://tinyurl.com/5yh6m6k5 • Movie S11: Video of rotation about 3D plot in manuscript Fig. 5(a) which shows the steady-states of experi-ments. The video starts from the perspective of the Fig. 5(a), then rotates about it to help readers understandthe 3D structure of the steady-states in the state-space of the experiments. It can be seen clearly that the Otherstates are distributed anisotropically about Trot, such that they are more widely spread across a plane spannedby the two linear, orange invariant manifolds, ( D s , https://tinyurl.com/2ykm5796 • Movie S12: Video of rotation about 3D plot in manuscript Fig. 5(b) which shows the steady-states of heteroge-neous simulations. The video starts from the perspective of the Fig. 5(b), then rotates about it to help readersunderstand the 3D structure of the steady-states in the state-space of the heterogeneous model. The Otherstates are distributed anisotropically about Trot, such that they are more widely spread across a plane spannedby the two linear, orange invariant manifolds, ( D s , https://tinyurl.com/3mvttf4d (b) View of experimental feature containing multiple networks (c) View of experimental feature with applied lightReentrant window Isolated 4-ring network Actinic inhibiting light(a) View of isolated4-ring network i) 4-ring networkii) 4-ring network with inhibited boundary and scalebar iii) 4-ring network depth Hiv) PDMS underreactors PDMSGlass d Ld m L m PDMSGlass
FIG. 1. (a) i) A single network of four coupled BZ reactors. The side lengths, L , of the reactors are 62 µ m, side-to-sidedistances, d , between reactors are 26 µ m.ii) The distance between a reactor and the inner surface of the inhibited boundary, d m , was 130 or 190 µ m. The thickness of the inhibited boundary, L m , was 100 µ m. iii) 50X magnification profilometer imageof cut microfluidic reactor network showing the out of plane height H of the reactors is 30 µ m and iv) 100X magnificationprofilometer image of cut microfluidic reactor network showing there is below 2 µ m of PDMS between the bottom of the reactorsand the glass below it. (b) View of running experiment with nine parallelized 4-node ring networks running simultaneously.The reentrant window is a piece of glass used to seal the BZ into the lattice. The black region along the borders, outside thereentrant window, is filled with oscillatory BZ reaction [1]. (c)
View of same experiment except with actinic, inhibiting lightselectively patterned onto the orders. . II. EXPLANATION OF BEST FIT MODELA. Fitting to experiments/choice of free parameters
The best fit model corresponds to substantially weaker than expected diffusive coupling. The diffusive couplingrate predicted by quasi-static diffusion between rectangular reactors was previously derived by Norton et al. [2]: k ideal = P DLd (1)where P is the partition coefficient, D the diffusion coefficient, L the side length of the reactors, and d is the distancebetween reactors [2, 3]. In our case D = 10 − [m s − ], P = 2 . L = 62[ µ m], and d = 26[ µ m] thus the idealizedcoupling rate is k ideal = 2 . [s − ]. We see the best fit value of k (manuscript Table V) is two orders of magnitudebelow k ideal . In past works with emulsion droplets the best fit value of k has been closer to one order of magnitudelower than expected [2–4]. Relative these past experiments with emulsion droplets in glass capillaries, densely packed (b) Sealing components: HeatsinkTECThermistor port Single feature ThumbscrewClamp frameO-ringElastic layerReentrant PDMS coated glass-slide with multi-ple networks (d) Experimental clamp:(c) Temperature control components:
TEC Clamp frame topClamp frame bottomGlass slideAluminium Aluminium PDMS layerO-ringReentrant Elastic supportGlass supportSingle sealed fea-ture of networks HeatsinkThermistor portThumbscrew Clamp frame bottomPDMS coated glass-slide with multiple networksClamp frame topThermal cutoff TECGlassPDMS (a) Sealing sets of BZ networks: . mm . mm μ m FIG. 2. Temperature controlled BZ clamp. (a)
CAD diagram of sealing of a set, or feature, of BZ networks. The networksare sealed above and below by glass. A volume of extraneous BZ is also sealed, however is irrelevant to the dynamics of thenetwork. Sealing is achieved by gradual adjustment of fine-threaded (6-80) thumbscrews. These apply force via a metal clampframe, and softened by an 2.5mm thick layer of PDMS. (b)
CAD diagram rendering of sealing components of the clamp. (c)
CAD diagram of the temperature controlling elements of the clamp. Temperature measurement occurs through a thermistor.The clamp is heated and cooled using 2 thermoelectric cooler (TEC) peltier elements placed symmetrically about the sample.Feedback between them is performed using an Arduino’s standard PID libraries. (d)
Photograph of temperature-controlledclamp containing a sealed network. All components of the clamp are clearly visible except the TECs, elastic layer, reentrantwindow, and o-ring. Unlike in simplified diagrams (a),(b), and (c), three other necessary components of the clamp are pictured:wiring, water circulating tubes connected to the heat sinks, and thermal cutoffs. The wiring connects the peltier TECs andthermistor to the controlling Arduino. The tubes circulate water through the heat sinks and cools the water down using fanradiators. The thermal cutoffs in series with the TECs short-circuit reversibly if the clamp become hot, above 50 ° C, by someerror. A meter stick is shown for reference (units of cm). emulsion droplets, and droplets in silicon chips, the experiments performed here had a much larger amount of oil-phasein the form of PDMS surrounding each 4 ring network Fig. 1. As Br is known to partition into and react with oiland PDMS, more PDMS would cause a reduction in inter-reactor coupling. We hypothesize this is the cause of thisdiscrepancy.Since the coupling rate equation Eqn.1 is composed of known parameters L and d , precisely prescribed by thephotolithographic process when the chip is fabricated, and the product of the two less certainly known quantities P and D , we can consider these best fits as a measurement of the product P D . Given the best fit k = 1 . × − [s − ], P D = 29 . µ m s − ].Further, the best fit model corresponds to only slightly heterogeneous reactors, naively making it a near ideallysymmetric network. To gain a sense for how heterogeneous the system is, we calculate the percent coefficient ofvariation in intrinsic frequencies of our best fit model. We consider a mean intrinsic frequency ω o = π [rad s − ].We can then consider the probability of an intrinsic frequency about the mean using the probability density func-tion(PDF) ρ (∆ ω ij ) = ρ ( ω i − ω o ) of our best fitted intrinsic frequency differences, manuscript Fig. 6(c). Noting thatthe distribution of intrinsic frequencies ρ ( ω i ) is then the same Laplacian PDF as ρ ( ω i − ω o ), just shifted to be about ω o , we can see the variance σ of intrinsic frequencies is the same as the variance of the PDF ρ (∆ ω ij ). We thenexplicitly calculate the variance of a Laplacian PDF, as it given by σ = 2 b where b is the rate parameter, so thestandard deviation is σ = √ b = √ × π × − rad s − (manuscript Table V). Using this and the mean value ω o , wecan compute a percent coefficient of variation 100 ∗ σω o = 10 √ × × − (3 × ) − = √ × . . LCD arraySampleGreen Kohler illuminationImage of LCD array x ob j e c t i v e i n f i n i t y - c o rr e c t ed - B ea m s p li tt e r T ube l en s CCD i) Kohler illumination arm : G r een L E D Len s A pe r t u r e 515n m f il t e r A pe r t u r e Len s Inverted stock projector lensDSLR infinity focused lensReduced image of LCD array (a) Ray diagram of programmable illumination microscope (b) Programmable illumation exampleProjected image (c) Picture of programmable illumination microscopeCCD view
LCD projectorDSLR infinity focused lensKohler illuminationTube lensCCDSampleii) Light perturbation arm: Inverted stock projector lens
FIG. 3. (a)
Ray diagram of a both i) microscopy imaging of BZ networks and ii) their spatiotemporal perturbation with light.i) Imaging of the networks functions through 515nm green transmission illumination, ii) Spatiotemporal perturbation is achievedby imaging a commercial liquid-crystal display (LCD) projector’s LCD array set to project pure blue light onto the sample.The blue light efficiently excites the photo-sensitive, oscillation-inhibiting catalyst. (b)
Example of how spatiotemporal lightpatterning is achieved and its quality. Images are projected onto the sample by sending images to a VGA/HDMI port connectedto the projector, an example is labeled ’Projected image’. This can be done through powerpoint, however we use MATLAB.This image then hits the sample, homogeneous and in crisp focus, shown in ’CCD view’. (c)
Picture of full experimentalapparatus, including a clamped, temperature-controlled sample. Overhead light shown in picture not on during experiments. derived results. How oscillators respond to heterogeneity doesn’t just depend on coupling strength k , it also dependson strength of the interaction function H . Letting H max be the maximum absolute value of H , in a pair of oscillatorsthere is a critical nondimensionalized intrinsic frequency | ∆ ω kH max | = ∆ ω ∈ [1 ,
2] at which no phase-locked steady-statecan exist past, initiating phase slipping. As the pair is symmetric when ∆ ω = 0, the initiation of phase-slippingat finite ∆ ω is a point where the system is far from the fully symmetric limit. Although we are studying a largernetwork, so ∆ ω > ω inmanuscript Fig. 6(c) allows us to estimate how far experiments are from being dynamically ideally symmetric. Wesee from this plot there are experiments which sample the tail of the exponential distribution of | ∆ ω ij kH max | and arefar from truly symmetric dynamics close to | ∆ ω ij kH max | = 1. Thus, although the reactors possess very similar intrinsicfrequencies, the interaction function H for our coupled BZ microreactors is sufficiently weak that experiments aredynamically not ideally symmetric. As a result, we can understand why experimental steady-states were at timesdistant from their ideal, H/K predicted values in manuscript Fig. 5(a) and in Fig. 6(a)(b), in spite of the reactors’small frequency dispersion.The necessity of excitatory coupling is due to observation of in-phase synchronization during experiments. Withoutexcitatory coupling, k e = 0, best fits to experiments with initial conditions that converged to Pronk and Bound wereinaccurate as they would always converge to Trot. Although with no excitatory coupling, k e = 0, the phase model hasPronk and Bound as attractors, the basins of attraction of Pronk and Bound are much smaller than experimentallyobserved. Allowing finite excitatory coupling, specifically k e = .
05, achieved excellent fits to experimental data.Coupling across the diagonal was observed to have no impact. Surprisingly, f at nearly any value from 0, a perfectring with no cross-talk, to 1, a fully connected network, did not have a strong impact on the error of the fitting.Investigating the Pronk, Bound, Rotary Gallop, and Trot states, we found the stabilities of these states are only veryweakly altered by the value of f . For this reason, f = 0 was selected. B. Best fit model:All steady-states
In addition to the analysis of the best fit model included in the manuscript, we attempted to find all its steady-states: ¯ θ such that ¯Ψ(¯ θ ) = ¯0. This was achieved through dividing state-space into a NxMxM grid and initializingNewton-Ralphson rootfinders on ¯Ψ(¯ θ ) with no heterogeneity, ∆ ω ij = 0, at the sites of that grid. The result ofmany iterations of the rootfinding ¯ θ n was considered a valid steady-state if ¯Ψ(¯ θ n ) . ¯Ψ(¯ θ n ) < × − . The identified S1 S2S1 S25 10 15 20 25Time[ks] P ha s e d i ff e r en c e [ r ad ] L i gh t[ A . U ]
543 10 11 12 P ha s e d i ff e r en c e [ r ad ] L i gh t[ A . U ] (a) Full experiment(b) Beginning of frequency locking (c) Beginning of phase-lockingPhase-lockedFrequencies similar1 234 FIG. 4. (a)
View of a single 4 reactor network over time during an entire experiment. The top plot is of light transmittedthrough the reactors. It is normalized from the raw intensity measured as an average of a square of pixels inside a reactor[Fig. 1(b)]. Specifically if the raw signal was I ( t ), with a maximum value during the experiment was I max and a minimumat I min , the normalized intensity was ( I ( t ) − I min ) / ( I max − I min ). As time goes on, the amplitude of the light peaks goesdownwards, while the baseline increases. The lower panel tracks the phase differences between the reactors versus time.Eventually the reactors reach a phase-locked steady-state in which all reactors are firing antiphase their two neighbors. Sincethe phase differences between the reactors are never perfectly stationary, we use a specific threshold for when phase differencesevolve slowly enough to be considered phase-locked in Appendix B. At about 20,000s the experiment ceases to be consideredphase-locked as its phase differences become noisy. (b) An inset of early times in the experiment. Note that the frequenciesand amplitudes of oscillation of the reactors are initially not the same. For example, the interval between the first two spikes ofreactors 2(red) and 4(green) are different, and the amplitudes of 3(blue) and 4(green) are different. However, at about 4,400sthe reactors possess similar frequencies by the threshold in Appendix A Protocol. (c)
An inset of later times in the experiment.Note that the frequencies and amplitudes of oscillation of the reactors are very similar. At the green line marking 12,050s, thethreshold for phase-locking is met. steady-states are shown as space-time plots with their Lyapunov exponent and topological index in Fig. 6 and ascolored points in the state-space in Fig. 7. The Lyapunov exponents were calculated as the maximum real eigenvalueof the semi-analytically computed the Jacobian of the dynamics [Eqn.2] at a given steady-state, utilizing numericallycomputed derivatives of the interaction function, H , manuscript Fig. 6(a). The state-space contains 168 steady-statesFig. 6; 10 of which are attractors or repellors while the remaining 158 are all saddle steady-states. We believe thismethod may have discovered all steady-states since all indices sum to 0, the Euler characteristic of the state-space.Regardless of whether we found all steady-states, it does show: (1) The unique stable steady-states are the H/K point G be s tf i t ( x ) [ r ad ] χ[rad] Best fit BZ reactor pair interaction
FIG. 5. The pairwise interaction function between two BZ reactors coupled by the best fit interaction function: ddt θ = kG ( θ ) = k [ H ( − θ ) − H ( θ )]. We can see that two negative slope zero crossings occur at θ = 0 , π - thus the best fit modelhas both in-phase and antiphase as stable attractors. invariant manifolds and (2) the state-space possess a wealth of unstable, saddle steady-states. This method was tooslow to be used to efficiently calculate the steady-states of the model with heterogeneity. III. COMPUTATION STABILITY OR TRANSVERSE LYAPUNOV EXPONENTS OF THEINVARIANT MANIFOLDSA. Example: computing transverse Lyapunov exponents of an invariant manifold
Following the algorithm present in Appendix E for the invariant manifold ( D p , D p ):1. We choose the invariant manifold ( D p , D p ). It is a 2D manifold in the 3D state-space of the 4 ring network.2. The invariant manifold ( D p , D p ) is parametrized by ¯ θ = ( f , − f , f ), manuscript Table II. The set of tangentvectors which span it are T = {√ − (1 , − , , (0 , , } . The one vector normal, which we put in our set ofpossible normal vectors N = { − (1 , , } .3. We now compute the unitary transformation matrix P . In the case of the 4 ring’s invariant manifold this isvery simple: P = [¯ t ¯ t ¯ n ] T = 1 √ − √
21 1 0
4. For our particular 4 ring network the Jacobian can be analytically computed: H ( x ) = dH ( χ ) dχ | χ = x , J (¯ θ ) = − k " [ H ( θ )+ H ( − θ )+ H ( θ + θ + θ ) ] [ − H ( θ )+ H ( θ + θ + θ ) ] H ( θ + θ + θ ) − H ( − θ ) [ H ( θ )+ H ( − θ ) ] − H ( θ ) H ( − ( θ + θ + θ )) [ H ( − ( θ + θ + θ )) − H ( − θ ) ] [ H ( θ )+ H ( − θ )+ H ( − ( θ + θ + θ )) ] + f k H ( θ + θ ) [ − H ( θ + θ ) + H ( θ + θ )] − H ( θ + θ ) H ( − ( θ + θ )) [ H ( − ( θ + θ )) + H ( θ + θ )] H ( θ + θ ) − H ( − ( θ + θ )) [ H ( − ( θ + θ )) − H ( − ( θ + θ ))] H ( − ( θ + θ )) (2)Note that since we numerically computed H , H is computed numerically as well. ����� - ���� ���� �� [ � / � ] ���� { θ �� � θ �� � θ �� � θ �� } ����� ������� - { } - - { } - - { } - - { } - { } { } { } - { } - { } - { } - { } - { } - { } - { } - { } - { } { } { } { } { } { } { } - { } - { } { } { } { } { } { } - FIG. 6. A list of numerically discovered symmetrically unique steady-states of the best fit model of the 4 ring network withoutheterogeneity. For each state we computed: (i) An example of the space-time plot of the state, (ii) The Lyapunov exponent ofthe steady-state with units [s − ], indicating the exponential rate at which nearby trajectory collapse to, if λ <
0, or divergefrom, if λ > π , i.e. as fractions of a period (v) the topological index of the steady-state and (vi) the numberof distinct states in state-space equivalent the unique state on the left. To be precise the number of distinct states is equal tothe number of 4-tuples of phase differences ( θ , θ , θ , θ ) generated by operating on the state on the left by every elementof the symmetry group D and then modulo 2 π . - - - Pronk RightRotary GallopBound Pace Trot 2-Saddle Pace-Trot separatrixPronk-Pace separatrixAttractorRepellor1-Saddle TransverseLyapunov exponent [1/ks]
FIG. 7. All steady-states of best fit model without heterogeneity which are in the half state-space θ > π − θ , colored bytype: Blue are Attractors with index=-1, Orange are saddles with 2 attracting direction saddles and thus index=1, Black aresaddles with 1 attracting direction and thus with index=-1, Red Repellors with index=1. Points are faint if they are a periodicinstances of another steady-state. The unstable steady-stats which form the separatrices between Pronk and Pace, and Paceand Trot are also marked.
5. We compute J = P.J.P − . We let f = 0 for simplicity. J = − H ( − θ ) − H ( θ ) − H ( θ ) − H ( θ ) − H ( θ ) − H ( − θ ) − H ( θ ) ( H ( θ ) − H ( − θ ))0 0 − H ( θ ) − H ( θ )
6. We extract the ’block’ of dynamics corresponding to perturbations in the only normal direction[ − H ( θ ) − H ( θ )]7. In general, this extracted block will not be 1x1 and the max Transverse Lyapunov exponent will be its maxreal eigenvalue. Since the matrix is 1x1 the value itself is the MTLE. We conclude that a infinitesimally smallperturbation normal the ( D p , D p ) invariant manifold δn will evolve according to the equation ddt δn = [ − H ( θ ( t )) − H ( θ ( t ))] δn .0 B. Maximum transverse lyapunov exponents of all weakly coupled 4 ring networks
We now present the MTLE for a 4 ring network of phase oscillators for a generic phase model with any interactionfunction H . With the exception of the linear invariant H/K manifold ( Z , J to identify the dynamics transverse to the manifold in J . The result of these calculations are shown in Table I is theblock of J corresponding to components normal or transverse the each invariant manifold. TABLE I. Decoupled transverse linear dynamics of invariant manifolds universal to four-rings( f = 0). H ( x ) = dH ( χ ) dχ | χ = x .Note that the max transverse Lyapunov exponent of an invariant manifold is the maximum real eigenvalue of the matrix in thecase of the linear invariant manifolds and a scalar in the case of the planar invariant manifoldsH K Transformed Jacobian J normal dynamics blockLinear invariant manifolds: D s D s − (cid:20) H (0) + H ( − θ ) H ( − θ ) H ( θ ) 2 H (0) + H ( θ ) (cid:21) D s − (cid:20) H ( π ) + H ( − θ ) H ( − θ ) H ( θ ) 2 H ( π ) + H ( θ (cid:21) Z (cid:20) J ( z , J ( z , J ( z , J ( z , (cid:21) Planar invariant manifolds: D p D p − H ( θ ) − H ( θ ) Terms for ( Z ,
1) invariant manifold: J ( z , = (2 H ( π + θ ) + H ( − π + θ ) − H ( θ ) − H ( − θ ) − H ( − θ + π )) − H ( − θ − π )) J ( z , = − ( − H ( π − θ ) − H ( − π + θ ) + H ( π − θ ) + 2 H ( − π − θ )) J ( z , = ( − H ( π + θ ) − H ( − π + θ ) + 3 H ( θ )) J ( z , = − (2 H ( π + θ ) + H ( − π + θ )) C. Transverse stability of invariant manifolds: Kuramoto model compared to BZ phase model
Having both parametrized guaranteed invariant manifolds (manuscript Table II) and provided firm expressions forthe MTLE of 4 ring networks’ invariant manifolds (Table I), we compare the invariant manifolds of 4 ring Kuramotomodel networks and those of the best fit model of coupled BZ reaction-diffusion network. In previous works somesimilarity was seen between the simplified Kuramoto model and inhibitory coupled BZ reactors in hexagonal lattices[5, 6]. The MTLE for networks of 4 nodes with perfect dihedral 4 symmetry (A) Kuramoto 4 ring, model with H ( χ ) = − sin( χ ) , f = 0 and k = 1, and (B) BZ best fit model [manuscript Table V] are shown as heatmaps appliedto H/K invariant manifolds in Fig. 8.We observe that the two 4 node networks are qualitatively different. The Kuramoto model’s invariant manifoldsare exclusively repelling except of ( D s ,
1) and ( D p , D p ) about Trot. The best fit BZ model’s invariant manifolds arealmost exclusively attracting, with the exception of ( Z , IV. EULER CHARACTERISTIC OF A 3-TORUS
The Euler characteristic of a solid Polyhedra is given by: χ = V − E + F − C which depends on integer number of vertices V , edges E , faces F , and number of 3D objects C . The 3-Torus is acartesian cube with periodic boundary conditions. Thus it is a cubic polyhedra, which very few unique vertices, edges,and faces.1 - - - - - - (a) Stability of Kuramoto network invariant manifolds (b) Stability of BZ network invariant manifolds Transverse Lyapunov exponent[1/ks] Transverse Lyapunov exponent[1/ks]
FIG. 8. (a)
Maximum transverse Lyapunov exponent of 4 ring of inhibitory Kuramoto oscillators invariant manifolds. The( D s , D s ) and ( Z ,
1) invariant manifolds are fully repelling. The ( D p , D p ) invariant manifolds are half repelling and halfattracting. The ( D s ,
1) manifold is largely attracting. (b)
Maximum transverse Lyapunov exponent of best fit model of 4 ringof inhibitory BZ oscillators. All invariant manifolds are largely attracting, except ( Z , Vertices: There are 8 vertices of a normal cube. If we consider any vertex at the edge of the periodic cube, we seethey are related to one another by periodic boundary conditions. V = 1.Edges: There are 12 edges of a normal cube. Each edge of the periodic cube is repeated 3 times. For example theedge ((0 , , , (0 , , π )) is equal to ((2 π, , , (2 π, , π )),((0 , , , (0 , π, π )), and ((2 π, π, , (2 π, π, π )). Thereare thus 3 unique edges E = 3Faces: There are 6 faces to a normal cube. In the periodic cube each face is repeated twice. F = 3 Thus the Eulercharacteristic is 0: χ = 1 − − [1] T. Litschel, M. M. Norton, V. Tserunyan, and S. Fraden, Lab on a Chip , 714 (2018), arXiv:arXiv:1711.00444.[2] M. M. Norton, N. Tompkins, B. Blanc, M. C. Cambria, J. Held, and S. Fraden, Physical Review Letters , 148301(2019).[3] N. Tompkins, N. Li, C. Girabawe, M. Heymann, G. B. Ermentrout, I. R. Epstein, and S. Fraden, Proceedings of theNational Academy of Sciences , 4397 (2014).[4] N. Li, J. Delgado, H. O. Gonz´alez-Ochoa, I. R. Epstein, and S. Fraden, Physical Chemistry Chemical Physics , 10965(2014).[5] N. Li, N. Tompkins, H. Gonzalez-Ochoa, and S. Fraden, European Physical Journal E , 1 (2015), arXiv:15334406.[6] M. Giver, Z. Jabeen, and B. Chakraborty, Physical Review E83