Decoding Network Structure in On-Chip Integrated Flow Cells with Synchronization of Electrochemical Oscillators

The analysis of network interactions among dynamical units and the impact of the coupling on self-organized structures is a challenging task with implications in many biological and engineered systems. We explore the coupling topology that arises through the potential drops in a flow channel in a lab-on-chip device that accommodates chemical reactions on electrode arrays. The networks are revealed by analysis of the synchronization patterns with the use of an oscillatory chemical reaction (nickel electrodissolution) and are further confirmed by direct decoding using phase model analysis. In dual electrode configuration, a variety coupling schemes, (uni- or bidirectional positive or negative) were identified depending on the relative placement of the reference and counter electrodes (e.g., placed at the same or the opposite ends of the flow channel). With three electrodes, the network consists of a superposition of a localized (upstream) and global (all-to-all) coupling. With six electrodes, the unique, position dependent coupling topology resulted spatially organized partial synchronization such that there was a synchrony gradient along the quasi-one-dimensional spatial coordinate. The networked, electrode potential (current) spike generating electrochemical reactions hold potential for construction of an in-situ information processing unit to be used in electrochemical devices in sensors and batteries.

Scientific RepoRts | 7:46027 | DOI: 10.1038/srep46027 electrical potential perturbation of one electrode results in electrogeneration of chemical species, which can diffuse to the other electrode for collection and possibly new electrical stimulation 19 . With different mechanism, bipolar electrochemistry 20,21 can be used in a single flow channel to fabricate simple logic gates with inputs as voltage sources and outputs as optical signal from electrogenerated chemiluminescence 22,23 .
In this paper, we explore the coupling topology that emerges through the potential drops in the flow channel in a commonly used electrochemical lab-on-chip device, a multielectrode detector system in a microfluidic flow channel. We employ a complex, electrochemical reaction (electrodissolution of nickel) 24 that generates oscillatory current (and electrode potential) at constant circuit potential. Analysis of the phase of the oscillatory reaction provides an effective means for characterization of the spiking patterns due to the coupling through network topologies. The network topology is decoded using dynamical measurements of reaction rates and phase model analysis (connectomics). The impact of cell geometry on widely different (positive and negative, uni and bidirectional) coupling topology is interpreted with a theoretical model of the potential drop in the flow channel. The effects of the unique, position dependent network topologies on the features of the spatially organized partially synchronized states are analyzed as a function of electrode number (two to six) and cell geometry (position of the electrodes).

Results
Our general strategy to explore the electrical coupling among the electrode is as follows. First, we consider cell geometries defined as number of working electrodes, and placement of reference/counter electrodes. (Figure 1 shows the schematics of the three considered configurations with two working electrodes). A typical device has a flow channel in which a number of electrodes are placed. The chemical reactions take place on the surfaces of the electrodes, and the rate of the reaction strongly depends on the local electrode potential drop that drives the reaction. The electrochemical reaction generates current, that flows to the counter electrode placed in the reservoir. A potentiostat sets the potentials of the electrodes, such that the potential differences between the working and reference electrodes are constant. Because the current flow between the electrodes introduces potential drops through the electrolyte, the reaction rates measured on the electrodes depend on each other, and the dependence is affected by the spacing of the electrodes and the placement of the counter electrode (which defines the direction of the current) and the reference electrode (which defines a potential with respect to the electrodes are polarized). Therefore, our next step is the application of equivalent circuit modeling to obtain a mathematical formula for the coupling between the electrodes as a function of the total resistance of the cell and geometry parameters (e.g., resistance, which is proportional to the distance, between the electrodes) 25,26 . (Note that because the objective is to explore the effect of electrical coupling between the electrodes, we have chosen a chemical reaction, the electrodissolutin of nickel in sulfuric acid 24 , which is under kinetic control and thus not affected by mass transport effects). Finally, experiments are presented with spiking patterns of the oscillatory reactions to demonstrate the effects of the coupling topology.
The dynamics of the electrode potential of a single electrode in the array can be described with the use of a Randles equivalent circuit. The current (I k ) generated by each electrode that has capacitance/surface area C d is obtained from double layer charging and Faradaic current: where A is electrode surface area, E k is electrode potential, and J F,k is the Faradayic current density of the k-th electrode. By rearranging equation (1), In the presence of many electrodes in the flow channel, the currents generated by the electrodes contribute to the ohmic (IR) potential drops in the flow channel; these ohmic drops affect the electrode potentials of the electrodes in a non-trivial manner and result in network interactions among reactions.
The differential equations for the dynamical evolution of the electrode potentials can be written in a general form as: where V is circuit potential, and K 2→1 and K 1→2 are the coupling strength from electrode 2 to 1 and electrode 1 to 2, respectively. The derivation of the equations for the considered cell geometries can be found as Supplementary Note. We consider the same circuit potential (V) and total cell resistance (R 0 ) for each electrode; therefore the system parameters for the surface reactions are identical. Such approximation greatly simplifies the extraction of interaction topology between the reaction units. We assumed that the potential drop in the flow channel is nearly linear and thus the the resistance in the flow channel is proportional to the channel length. In other words, we assume that Ohm's Law is obeyed for the potential drop in the cell. This assumption was validated in numerical simulations and experimental measurements of the resistance of the flow channel 25,27 . Note that such assumption is valid when the quasi-one-dimensinal flow channels are much longer than the size of the electrode; the approximation is usually acceptable in a typical flow channel length of few mm and electrode sizes in the few μm range.
The coupling strengths obtained from equivalent circuit analysis between two working electrodes are given in Fig. 2. Ipsilateral (traditional) placement configuration. First we consider two working electrodes in the flow channel at ipsilateral (traditional) placement in which the reference and counter electrodes are at the same downstream end of the flow channel as shown in Fig. 1(a).
The theoretical analysis (see Fig. 2) showed that . The dynamics is affected by bidirectional positive coupling that depends on the total cell resistance R 0 , distance from downstream electrode to the reservoir (through solution collective resistance R C ), and electrode surface area A. Thus, the mathematical analysis suggest the existence of bidirectional positive coupling within the cell.
We performed experiments in a dual electrode setup to confirm the theoretical findings. In this configuration, the oscillatory metal dissolution of nickel in sulfuric acid exhibited in-phase synchrony (see Fig. 3(a)).
With an oscillatory system close to the onset of oscillations, the synchronization pattern is typically either in in-phase or anti-phase configuration. Positive coupling implies that electrode potential difference will create a cross-current between the electrodes, which diminishes the potential difference between the electrodes. Therefore, very often positive coupling between the electrodes will result in in-phase configuration. Indeed, in previous investigations with macrocells, where positive coupling was induced with an external resistance interface, in-phase synchrony behavior was observed 28 . In contrast, when negative coupling is present (e.g., with IR compensation), anti-phase oscillations can occur 28 . The in-phase synchronized oscillations (shown in Fig. 3 ) coupling topologies.
Extracting coupling topology from a network of dynamical units is a challenging task that requires advanced signal processing techniques 29,30 . Here we apply a method that uses phase models to characterize the interaction topology 31,32 . For uncoupled elements, the phase of the oscillations (φ) increases linearly with time at a rate of the natural frequency (ω), i.e., dφ/dt = ω. However, when the oscillations are coupled, there will be deviations from the linear variations: the deviations depend on the phase difference between the coupled elements, therefore, we can express the rate of phase change with the following equations 31 : where t is time, φ k and ω k (k = 1, 2) are phases and natural frequencies of the oscillations, Δ φ is the phase difference, and H k are the phase interaction functions. The interaction functions H k can be determined by plotting dφ k /dt − ω k vs. Δ φ from an experiment with phase slipping behavior 32 . As shown in Fig. 3(d), the interaction functions for for both units are nearly sinusoidal functions. Coupling topology extraction is based on the assumption that the coupling strength is proportional to the amplitude of the interaction function. Therefore, when the amplitude of the interactions functions of two dynamically identical units are the same, the coupling is symmetrical. Conversely, with only unidirectional coupling the phase dynamics of the driven node is affected, i.e., the amplitude of the interaction function of the driver node must be zero. The symmetrical and unidirectional couplings represent two extremes by which two units can interact. The coupling between two units can be described by the coupling strengths in both directions (i.e., K 1→2 and K 2→1 ), or by a combination of a mean coupling strength and a directionality index, We can obtain the directionality index (d) from the comparison of the amplitudes of the (first) Fourier harmonics (c k ) of the interaction functions H k : 1 . d = ± 1 and 0 indicate strictly unidirectional and bidirectional coupling, respectively. The interaction functions for the two electrodes at ipsilateral (traditional) placement ( Fig. 3(d)) show approximately positive sinusoidal functions with amplitudes c 1 = 0.05 rad/s, c 2 = 0.04 rad/s (for electrodes 1 and 2, respectively), therefore, the analysis supports the presence of positive bidirectional (d = 0.1) coupling between the electrodes (shown in Fig. 3(g)). Contralateral placement configuration. In another configuration with contralateral placement of reference/counter electrodes, the reference and counter electrodes are placed at the opposite end of the flow channel where the two working electrodes reside (see Fig. 1(b)).
The corresponding equivalent circuit analysis (Fig. 2) shows that there is no coupling from the downstream to the upstream electrode and the coupling from the upstream to the downstream electrode is a negative value that depends on the resistance (and thus the distance) between two working electrodes R 12 , i.e., 12 . Experiments with contralateral placement configuration were carried out with a dual electrode system, where a 2 mm diameter Ni electrode placed upstream to the working electrodes served as a (pseudo) reference electrode. In this configuration the observed anti-phase oscillations (Fig. 3(b)) indicate the presence of negative coupling. The phase interaction function (Fig. 3(e)) is a negative sinusoidal function with amplitude c 2 = 0.23 for the downstream electrode and a nearly straight line with negligible amplitude c 1 = 0.01 for the upstream electrode. Consequently, unidirectional coupling (schematically shown in Fig. 3(h)) is quantified by the large value of the directionality index d = 0.9. Note that by a simple movement of the reference electrode from the downstream to the upstream positions changed the sign and directionality of the coupling in the network.
Dual-reference electrode configuration. In a third arrangement (dual-reference electrode system, see Fig. 1(c)), the upstream 2 mm diameter nickel electrode serves as the reference electrode for upstream working electrode (WE1), the downstream Ag/AgCl/3 M NaCl reference electrode serves as the reference point for downstream working electrode (WE2), and both of the working electrodes share the same 0.5 mm thick Pt wire counter electrode at the downstream end of the flow channel.
In this configuration, the equivalent circuit analysis indicated that there is a unidirectional positive coupling that depends on the distance from downstream electrode to the reservoir, i.e., K 2→1 = 0 and . Thus by changing the position and number of reference electrodes, the coupling topology can be tuned between bidirectional positive, unidirectional negative, and unidirectional positive. With this setup, in-phase oscillations (Fig. 3(c)) were observed, which is an indication of positive coupling. Although the in-phase synchronization would imply some resemblance to the first (traditional) configuration with bidirectional positive coupling, the phase interaction function analysis (see Fig. 3(f)) indicates that only the downstream electrode is affected by coupling (c 1 = 0.00, c 2 = 0.10), and thus the coupling is unidirectional (d = 1.0), as shown in Fig. 3

(i).
Switched dual-reference electrode configuration. Alternative cell geometry was studied with the dual-reference electrode system, but the corresponding reference electrodes of individual working electrodes were exchanged (switched dual-reference electrode configuration, Fig. 4(a)). In this configuration, the upstream nickel electrode serves as the reference electrode for downstream working electrode (WE2), the downstream Ag/AgCl/3 M NaCl reference electrode serves as the reference electrode for upstream working electrode (WE1).
There is a positive coupling in one direction (downstream to upstream working electrode) and a negative coupling in the other direction, i.e., Fig. 2). With a simple change of the position of the downstream electrode the dominating positive ( Fig. 4(a)) or negative (Fig. 4(c)) coupling can cause in-phase ( Fig. 4(d,h)) or out-of-phase ( Fig. 4(g,k)) synchrony, respectively. In configurations where the positive and negative couplings had comparable magnitudes (Fig. 4(b)) in-phase synchrony was difficult to achieve: in some experiments complex oscillations with 4:5 entrainment behavior ( Fig. 4(e,i)), and in other experiments slightly out of phase synchronized oscillations were observed (Fig. 4(f,j)). The experiments thus show a potential for creating a largely different coupling directionality and strength by a simple change of the working electrode position WE2.
Configuration with three and six working electrodes. With the use of more than two electrodes, a network of reaction units could be obtained. Here we consider only traditional placements of reference and counter electrodes placed at the end of the flow channel as shown in the schematics in Fig. 5(a). Based on theoretical circuit analysis (see Supporting Note), the network topolog y consists of global coupling between each electrode pairs and an additional local coupling 0 0 23 between the two upstream electrodes, where R C is the solution collective resistance (depending on distance from electrode 3 to the reservoir), and R 23 is the solution resistance between working electrodes 2 and electrode 3. The definitions of K global and K local show that the coupling strengths between each electrode can be tuned by changing R 23 and R C .
The experiments were carried out with three micro-wires as shown in Fig. 5(a). Here we consider intermediate coupling strength at which the current oscillations of two out of three electrodes are fully synchronized, while third electrode is in the phase slipping state. To characterize the extent of synchrony between a pair of oscillators we use an order parameter r based on phase difference 33 defined as When the distance between WE2 and WE3 (D 23 ) is smaller than the distance between WE3 and the reservoir (D 3R ), the global coupling will drive the synchrony; the coupling topology is shown in the top panel of Fig. 5. Theoretical expectation of globally coupled population of oscillators in the partially synchronized state is that the identity of the synchronized oscillators will be uniquely determined by their natural frequencies 34 , thus independent of their geometrical positions. We performed a series of experiments with different electrodes when D 23 < D 3R ; representative experiments are shown in Fig. 5(a-c). As shown in the synchronization matrices in Fig. 5(d-f), when the oscillators exhibited partial synchrony, the two synchronized oscillators could occur from any possible pairs (1-2, 1-3, or 2-3), thus the synchrony pattern does not depend on the electrode placement. These results are in agreement with global coupling topology: partial synchrony does not exhibit spatial order.
When the distance between WE2 and WE3 (D 23 ) is larger than the distance between WE3 and reservoir (D 3R ), the local (upstream) coupling between WE1 and WE2 will drive the synchrony; the coupling topology is shown in the top panel of Fig. 6. Therefore, it is expected that the two upstream electrodes will be more synchronized than the other electrode pairs. To confirm this expectation, we repeated several experiments when D 23 > D 3R (see Fig. 6(a-c)). The synchronization matrices in Fig. 6(d-f) always showed that the two upstream electrodes (1-2) are more synchronized than the other electrode pairs. We refer to this state as spatially organized partial synchrony (SOPS).
As an extension of the spatial pattern shown before, a set of electrodes in the flow channel would be expected to have strong coupling among upstream electrodes and less coupling among downstream electrodes. In Fig. 7, we show the results obtained with six approximately equally spaced electrodes in the partially synchronized state. The three upstream electrodes (electrodes 1, 2 and 3) are synchronized with same frequency (0.476 Hz), the downstream electrodes are less synchronized (electrodes 4 and 5 with frequency of 0.520 Hz and 0.534 Hz, respectively), and electrode 6 is not synchronized with even larger frequency (0.601 Hz) as shown in Fig. 7. Thus, the unique, position dependent network topology imposes a dynamical behavior with a gradient synchrony; such SOPS behavior seems to be the general feature of the linear array of microelectrochemical oscillator system.

Discussion
We demonstrated that the electrical coupling between the electrodes can generate networks with a variety of different topologies in microfluidic devices. Chemical connectomics, the decoding of network topology from dynamical measurements, was proven to be a useful tool to verify the predicted connectivities in two-unit systems. (For larger systems, advanced signal processing could be applied that considers indirect interactions 35 ). The mathematical formulas obtained for the coupling strengths hold promise for estimating the electrical cross-coupling in microfluidic devices in practical applications. In electroanalytical applications using traditional configurations, where the reference and counter electrodes are placed in a reservoir at the end of the flow channel, the distance between the most downstream electrode to the reservoir should be minimized to avoid global coupling. Attempts to eliminate coupling by moving the reference electrode to upstream position certainly facilitates the control of electrode potential of the upstream electrode, however, this configuration introduces dangerous negative (unidirectional) coupling between the electrodes. When this negative coupling is strong, it could create highly nonuniform reaction rates on the electrode surface, in particular, with reactions with negative slopes in the current vs. potential characteristics 36 . Therefore, when electrode arrays are employed, e.g., with scanning electrochemical microscopy 37 , care must be taken to allow current to disperse in the cell (e.g., by parallel placement of counter and working electrodes) to eliminate electrical coupling. The electrical coupling can be most adverse with generator-collector type of devices in linear flow channels 11 ; in thise case, small electrode areas and addition of parallel resistors could minimize the coupling effects.
The micro-electrochemical system studied here can accommodate large number of electrodes, and thus could represent an important experimental system for large set of networked chemical reactions. With the use of multichannel potentiostats, multiple groups of working electrodes could form modular networks, and the inter and intra-group coupling could be controlled by the application of multiple reference electrodes assigned to each groups. Such complex networks could be designed for specific applications. For example, a network with positive and negative connections can be designed with oscillatory units for associate memory arrays that can be used for pattern recogniation and storage 38 . By application of such information processing, the device could be integrated with traditional microfluidic sensors and batteries to create an all-electrochemical intelligent sensor device.

Methods
The experimental system is a microfluidic flow reactor in which 100 μm diameter nickel electrodes (embedded in epoxy) are sealed with polydimethylsiloxane chip. 2 M sulfuric acid/0.01 M NiSO 4 electrolyte is pumped in a 200 μm (width) × 100 μm (height) flow channel at flow rate Q = 1.5 μL/min that corresponds to a linear velocity of 0.125 cm/s. The Ag/AgCl/3 M NaCl reference electrode (exhibiting 209 mV electrode potential vs. standard hydrogen electrode) and a 0.5 mm thick Pt wire counter electrode are placed into the center of the reservoir. Further details about the the construction of the cell are given in a previous publications 25,26 .
When a single electrode is present in the flow channel and the electrode is polarized at a constant potential 1.6 V (vs. the reference electrode) the electrode reactions generate oscillatory dissolution rate (measured as the current) due to the secondary passivation of nickel 25,26 . The oscillation occur because of the interactions of a fast positive feedback loop and a delayed negative feedback loop in the reaction mechanism. Because of the negative differential resistance (decreasing current with increasing electrode potential), when the electrode potential increases, the current, and thus the ohmic (IR) potential drop in the cell decreases; as the circuit potential is kept constant, the electrode potential will thus increase which results in further decrease of the current. This fast positive feedback loop is then followed by delayed response of relaxation of oxide film coverage to steady state value; because the relaxation is relatively slow process compared to the time-scale of the variation of electrode potential, the feedback loops generate oscillatory metal dissolution. (We note that although we investigate the process in the oscillatory region, the derived networks are inherent properties of the cell configurations; we chose oscillatory reactions because with these reactions the dynamical behavior is is strikingly different for the different networks).
With more than one electrode, each of the electrodes exhibited smooth current oscillations (nearly sinusoidal waveform close to the Hopf bifurcation point) obtained at constant circuit potentials with 30 mV above Hopf bifurcation. (Because the onset of oscillations occurs at slightly different potentials for each experiment, we report the circuit potentials for every experiment).
The schematics of electrochemical cells are shown in Fig. 1(a-c). Ipsilateral (traditional) placement of reference/counter electrodes ( Fig. 1(a)) refers to the reference electrode and counter electrodes positioned at the same downstream end of the flow channel. For contralateral placement of reference/counter electrodes, the reference electrode is placed at the upstream end of the flow channel ( Fig. 1(b)). A 2 mm diameter Ni electrode, placed upstream to the working electrodes, serves as a reference electrode and a 0.5 mm thick Pt wire serves as counter electrode. Dual-reference electrode system is shown in Fig. 1(c), where the upstream 2 mm diameter nickel electrode serves as the reference electrode for upstream working electrode (WE1), the downstream Ag/AgCl/3 M NaCl reference electrode serves as the reference point for downstream working electrode (WE2), and both of the working electrodes share the same 0.5 mm thick Pt wire counter electrode. The switched dual-reference electrode configuration refers to cell that the upstream nickel electrode serves as the reference electrode for downstream working electrode (WE2), the downstream Ag/AgCl/3 M NaCl reference electrode serves as the reference electrode for upstream working electrode (WE1). Note that once the electrodes are integrated in the flow channel, it is very easy to switch between the different configurations by connecting the corresponding working, reference, and counter electrodes to the potentiostat.