Simulating waves, chaos and synchronization with a microcontroller

Simulations of pattern formation,^9–11 excitable media,^12–14 and oscillating systems in general require the solution of multiple coupled ordinary differential equations (ODEs) in time and partial differential equations (PDEs) in space.¹⁵ These systems can span from simple models such as FHN^6,7 and Hodgkin-Huxley (HH),¹⁶ consisting of two and four variables, respectively, to complex models with more than 80 variables.¹⁷ In many cases, the time and space discretization and the number of equations necessary to describe the system require the use of substantial computational resources such as supercomputers.^18–20 As computational power continues to follow Moore’s Law, which describes a doubling in computational speed about every two years, even 50 years after its conception,²¹ it has become easier to simulate large complex systems using common desktop computers.^22,23 In recent years, this computational gain has been extended with the use of graphics processing units (GPUs) as accelerators for scientific computing.^24–26 However, the increase in CPU power has also allowed the growth of another kind of resource to perform operations: microcontrollers (MCs).

The first MC that combined read/write memory, processor, and clock on one chip was created in 1971 by Texas Instruments. When introduced to the electronics industry in the form of the TMS-1000, it was used widely in calculators, toys, and games and sold at an affordable $2 per unit in bulk orders.²⁷ In response, Intel developed a MC for control applications in 1977, combining RAM and read-only memory (ROM) on a single chip. These two influenced the development of a range of applications, including personal computer (PC) keyboards. In the 1990s, advanced MCs were developed with electrically erasable and programmable ROM memories, which are still in use today by Atmel and Microchip. Today, MCs are even smaller, more powerful, and cheaper and are used in both active research and everyday consumer items such as phones, automobiles, and household appliances.^28,29

In a recent example, Serna and Joshi simulated the logistic map using LEDs driven by a MC.¹ They divided the $[0,1]$ interval domain of the logistic map into ten equal parts and mapped each to an LED that was powered when the trajectory of an initial seed lay in the corresponding interval. Fixed points were indicated by an LED that remained lit, while a group of such LEDs indicated a periodic orbit.¹ Separately, Mahmud et al.³ modeled 1D cardiac tissue propagation using analog circuits and digital signal controllers such as dsPIC MCs, wherein each cell’s voltage response was dictated by the Luo-Rudy phase I model.³⁰ MCs also have been used to simulate the HH neuron model with output via the serial port or to an LCD;³¹ in addition, field programmable gate arrays (FPGAs) have been used to solve the HH model.^32,33

Because of their flexibility and low cost, MCs are becoming an important tool for physics education in the class room.^34–36 For example, Soriano et al.³⁷ discuss the use of different MCs in robotics for classroom activities and in the advancement of automatic control.³⁷ Also, different courses such as Cornell University’s ECE 5760 “Advanced Microcontroller Design and System-on-Chip” have taught how to use a variety of MCs as components in electronic design. These lectures are listed on their website and go through multiple topics using a DE1-System-on-Chip (SoC) including visualizing the Mandelbrot set and simulating the Lorenz system. This course teaches more in-depth how the SoC works with other devices such as FPGAs and digital differential analyzers.³⁸ 

In this paper, we further demonstrate the use of MCs to simulate three biophysical systems known for their rich dynamics including stable and unstable traveling waves, synchronization, spiral waves, and chaos. We show the dynamics of waves in rings of neuronal tissue whose dynamics are governed by the FHN model.^6,7 The FHN is a simple model that is used as a semiquantitative description of the dynamical behavior of an excitable neuron and is similar to the Hodgkin-Huxley model.¹⁶ This model, which has become a central example in reaction-diffusion equations, uses only two variables: one that mimics the membrane potential and another that mimics a recovery variable.

We then use the 3-variable FK phenomenological model of cardiac dynamics.⁸ The three dynamic variables of the model are membrane potential and two gate variables; together, they govern the behavior of currents that represent the flow of $Na+$⁠, $Ca2+$⁠, and $K+$ ions across the membrane. This model produces a more realistic cardiac action potential (AP), i.e., the spike and dip of the membrane potential in cardiac cells that occurs when an electrical pulse passes through the tissue, than the FHN model. This model is capable of demonstrating complex oscillatory pulse dynamics known as alternans³⁹ comparable to those observed in cardiac experiments.^40–43 

Next, we use MCs to illustrate synchronization,⁴⁴ a spontaneous emergence of order that can be observed in coupled oscillators such as some cardiac cells⁴⁵ and fireflies.⁴⁶ We use the Kuramoto model,⁴⁷ which describes the phase advance of coupled oscillators and can reproduce synchronization of weakly coupled homogeneous oscillators. While useful in some systems, such as the collective action of cardiac cells, synchronization is also responsible for destructive phenomena such as epileptic seizures⁴⁸ and the resonance caused by people crossing the London Millennium Footbridge while inadvertently synchronizing their footsteps with the oscillation of the structure.⁴⁹ 

This paper is organized as follows: Sec. II discusses the setup of the Arduino microcontroller. Next, Sec. III discusses different implementations of the FHN including a single oscillator that can be excited via touchscreen, a 1D ring of oscillators producing waves in a 60-LED ring, and a 2D grid of oscillators capable of producing spiral waves through interaction with the LCD touchscreen. Section IV discusses the results of visualizing more complex waves produced in space by a period-two bifurcation using the FK cardiac cell model and the 60-LED ring. Then, Sec. V demonstrates the results of simulating the KM with global and nearest-neighbor coupling using a 64-LED grid as physical space and a 24-LED ring, and with nonlocal coupling using the 60-LED setup, where chimeras are exhibited. In Sec. VI, we show how a simulation on a microcontroller can interact in real time with a web browser. Finally, the potential advantages and improvements on the system are discussed. All the programs for the Arduino and the Java applets are given in the supplementary material.

Microcontrollers come in all sizes and architectures, some even with their own Linux kernel installed and ready to use.⁵⁰ For simplicity and functionality, we use an Arduino Due, which is based on a 32-bit ARM core MC (cost: $39.95). This unit is connected to the following different examples: (i) one Adafruit NeoPixel NeoMatrix square 8 $×$ 8 RGB LED array ($34.95); (ii) four NeoPixel 1/4 60 RGB LEDs ($9.95 for each part); (iii) one circle Adafruit NeoPixel 24 RGB LED ring ($19.95); (iv) one Adafruit $3.5″$ TFT 320 $×$ 480 + Touchscreen Breakout Board w/MicroSD Socket ($39.95), which also needs a breakout board for some more complicated connections; and (v) a photoresistor. The 60-LED ring comes in four sectors and must first be connected by soldering the parts together. One of these setups is shown in Fig. 1, and others are shown in the supplementary material. Other displays can be used if a touchscreen is not needed or if a smaller size, such as the $2.8″$ TfT LCD with Touchscreen Breakout board w/MicroSD Socket ($29.95), is desired.

The Arduino can be programed with a language whose syntax is similar to C and Java. The codes (called “sketches”) can be uploaded from a computer through a microUSB to USB cord using an Integrated Development Environment (IDE) that can run and control the LEDs even when not connected to the computer.

For detailed explanation of the Arduino connections for the following examples, see the supplementary material.

The FHN model is one of the simplest models used to describe excitable systems such as chemical reactions,⁵¹ neurons,⁵² and cardiomyocytes.⁵³ It uses two variables to govern its dynamics, one fast denoted as $v$ and one slow denoted as $w$⁠. Depending on the system described (e.g., chemical, cardiac-neural, and population dynamics), $v$ can be associated with the activator, voltage, or prey, respectively, and $w$ with the inhibitor, ionic current, or predator. The model is thus described by two coupled ODEs,

Parameters $a$⁠, $b$⁠, $d$⁠, $δ$⁠, and $ϵ$ affect the location and type of equilibria and the speed of oscillations.⁶ Depending on model parameters, dynamics can be oscillatory or excitable. In these reaction type models, a difference in time scale between the two variables $v$ and $w$ is achieved by using $ϵ≪1$⁠, in this case making the dynamics of $v$ faster than $w$⁠.

As a first example, we use the sketch “FHN_AP.ino” (see the supplementary material) to show the dynamics of a single FHN cell in an oscillatory regime. The $v$ variable oscillates and its value is plotted in time on an interactive LCD (Fig. 2). The LCD displays the current value of $v$ at the left side of the screen and propagates the value in time to the right with the right side being the furthest back in time.

If untouched, the screen will trace the auto-oscillatory signal, or action potential, in time as shown in Figs. 2(a) and 2(b). The FHN action potential is characterized by a fast upstroke, a plateau, and a fast down-stroke to the rest state. In this regime, the oscillator can still be perturbed and excited by a stimulus during its rest state.

To incorporate this functionality into the simulation, the value of $v$ can be made to increase by an amount proportional to the pressure p.z applied anywhere on the LCD touchscreen by simply adding the following line:

v += p.z/130,

where the maximum pressure the Arduino can recognize (1000) is scaled so that it corresponds to a value that is comparable to values that the oscillator can achieve normally. This code adds some value between 0 (no press) and $7.7$ (the hardest press) to the $v$ variable during the time of the interaction. However, this implementation can be modified depending on how the user want the pressure to relate to the cell excitation. This way, by touching anywhere on the screen, the cell will get excited not only proportional to the pressure but also for the duration of the user’s touch on the screen. Accounting for this interactivity modifies the FHN equations as follows:

where $I$ is the impulse the user gives the oscillator through pressing the screen.

Depending on when during the oscillation's period the screen is touched, the cell may excite if the stimulus is strong enough.⁶ If the cell is excited during the rest period, it can easily produce early activations, as shown in Figs. 2(c) and 2(d), where the excitation is shown on the left [Fig. 2(c)]. After the simulation continues, there is a small action potential as the oscillator then continues back to rest [Fig. 2(d)]. However, if attempting to excite the cell while it is already excited or in the refractory period, the oscillator will not excite completely [Figs. 2(e) and 2(f)]. Figure 2(e) shows a user interacting with the screen exciting the oscillator a few times. Figure 2(f) shows the results of the user’s additional excitations. The first produces another early activation, but for the other three, the oscillator is not able to excite from the refractory period resulting in small excitations that die out.

The complete sketch for this example is listed in the supplementary material as “FHN_AP.ino.”

In systems of coupled cells, pulse solutions, multiple pulse solutions, and periodic traveling waves in the FHN model have been proven to exist⁵⁴ and their stability has also been studied.^55,56 Specific studies of rotating wave solutions on a circular ring have shown that the existence of these solutions and their possible wavelengths can depend on the diffusion coefficients and the domain size.⁵⁷ This study focuses on these traveling waves, one of the main known solutions. They can be obtained by exciting a group of cells to a higher voltage potential and integrating forward in time.

When cells in a ring are excited, they will elicit waves propagating in each direction, i.e., clockwise and counterclockwise, from the stimulus site; however, by applying unidirectional block,⁵⁸ it is possible to start waves in only one direction. These rotating waves have been shown to occur in cardiac tissue and act as anatomically induced rotating waves, which are a type of reentrant wave, causing tachycardia.⁵⁹ In the FHN model, there is a small parameter regime where waves can be elicited intermittently in each direction⁶⁰ from a single stimulus, a state that can be considered a one-dimensional spiral wave.

Here, the FHN model [Eq. (2)] is modified to have spatial coupling by diffusion,

The parameter values $a=0.2$⁠, $b=0.5$⁠, $d=1$⁠, $δ=0$⁠, and $ϵ=0.0095$ are used to simulate a traveling wave. For this case, an explicit Euler integration method is used with values $dx2=1$ and $dt=0.4$⁠, which leads to convergent results in 1D.⁶¹ The solution is plotted in real time using a ring of LEDs, where each LED represents a cell in the simulation. Figure 3(a) shows propagation of a stable wave (red LEDs) along the LED ring with a period of a complete rotation of 3.27 s. Also, this wave has a wavelength of about 18 LEDs, or an arclength of $108°$⁠, depending on what brightness is considered as part of the excitation. These measurements are easy to take: a stopwatch or a video taken via a camera or a phone can be used to track the wave completing a cycle. For more accurate measurements, a tracking program can be used to monitor the wave front and back propagating around the LED ring. The Arduino Serial can also be used to print out the values of each cell.

For comparison, a simulation of the same system using a Java applet is plotted on the screen. Figure 3(b) shows the variable $v$ (white) in space at different times. For the same parameters as used in the Arduino, the wave speed and wavelength can be measured. For the wavelength, a threshold is set for the voltage and the number of cells that are above that threshold at a certain time is measured. Therefore, the wavelength will change slightly depending on the threshold.⁶¹ For easy calculations, a threshold is picked to be $v=0$ and a wavelength of 20 cells, or $120°$ for the arclength, is measured. One complete rotation takes 154 unit time steps.

Figure 3 compares the output for both the LEDs and the Java applet and also plots the values of the calculations [Fig. 3(c)], which can be read out using the Arduino Serial while connected to a computer. In both cases, it can be seen that once the wave is stable, the wavelength and amplitude are also stable.

The complete sketch for this example is listed in the supplementary material as “FHN.ino.”

It is also possible to connect a photocell to the Adafruit ring to have a simple interaction via incoming light. Doing so makes the simulation interactive by touching the photocell, and a pulse with voltage dependent on the amount of light blocked can be programed to propagate out in both directions along the ring. This may cause the wave that is traveling along the LED ring to break up. Connections can be made following the user’s guide by LadyAda on the Adafruit website.⁶² 

Similar to the change made for interacting with the FHN single oscillator with the touchscreen in Eq. (2), an impulse is added to one oscillator in the ring. In this case, the impulse is given as follows:

I = (1023 - photocellReading)/130,

where photocellReading takes a value in $[0,1023]$⁠. As the sensor gets darker, this line will make the impulse higher with the highest value about $7.9$⁠. Similar to the single cell, if the impulse is not large enough, the pulse will die, which will be seen by an LED briefly lighting up and then dimming.

Two examples of the photocell interacting with the FHN wave are shown to illustrate different dynamics. Figures 4(a)–4(c) show a user covering the photocell, which causes a cell on the right hand side of the 60-LED ring to excite. This excitation then produces a propagating wave. The wave can propagate in both directions, but here, the cells in the clockwise direction are still in the refractory period, so that it can only propagate counterclockwise. Because a wave was already initialized and traveling counterclockwise, the new wave collides with the original wave [Fig. 4(d)] and they annihilate each other. If left with no more interaction, there will be no more waves propagating around the ring [Fig. 4(f)].

In the second example, all cells begin at rest, shown by all the LEDs being off (Fig. 5). The user then covers the photoresistor, which excites the first cell [Fig. 5(a)]. Because none of the cells are currently in the refractory period, the wave propagates in both directions until the two wave fronts meet and collide [Figs. 5(b)–5(d)] and the wave annihilates itself [Figs. 5(e) and 5(f)].

The complete sketch for this example is listed in the supplementary material as “FHN_Photocell.ino.”

Reaction-diffusion systems display several emergent behaviors as the dimension is increased. In 2D, besides plane and circular waves, which can be seen as an extension of the 1D dynamics, more complicated dynamics like spiral waves can be formed when the symmetry is broken.⁶³ To illustrate the induction of spiral waves, we extend the FHN to a two-dimensional grid-lattice topology with the nearest-neighbor coupling. Spiral waves can be seen in many physical and biological systems, such as cardiac ventricular fibrillation⁶⁴ and retinal spreading depression.⁶⁵ The interaction of external agents with the spiral waves can also create more complicated dynamics.

For visualization and interactivity, the Adafruit $3.5″$ TFT 320 $×$ 480+Touchscreen Breakout Board w/MicroSD Socket is used with the Arduino. This device has $320×480$ pixels and is interactive via pressure. The Arduino Due is limited in memory, however, so the largest system simulated for this example is $64×96$ with each simulated cell being $5×5$ pixels in size. Larger domains can be implemented with more expensive MCs such as Raspberry Pi.

For this simulation, Eq. (3) is extended to 2D where $vi$ becomes $vi,j$ and $wi$ becomes $wi,j$⁠. Also, the Laplacian of cell $vi,j$ is now calculated by using its four nearest-neighbors.⁶¹ To ensure convergent solutions, $dt=0.06$ and $dx2=1.2$ are used with the Euler integration. Also, to generate rotating spiral waves, parameter values of $a=2.2$⁠, $b=1.0$⁠, $d=0.05$⁠, $ϵ=1$⁠, and $δ=0$ are used.

To simulate spiral waves, a particular set of initial conditions is used.⁶³ By setting the potential $v$ of half a row of cells and the recovery variable $w$ of an adjacent row on one side higher than the resting potential, conduction block, a term borrowed from cardiac dynamics, occurs, causing the wave to propagate in only one direction.

The pressure sensor of the Adafruit display is used again to make the simulation interactive. Touching the screen will activate the section of cells touched, increasing the potential and changing the local dynamics. The Arduino can determine the location of the touch and transfer that to a location in the simulation grid. A radius of interaction can be set. For the simulation to activate the area of cells, the following lines are included:

if (p.z > MINPRESSURE && p.z < MAXPRESSURE) {

p.x = map(p.x, TS_MINX, TS_MAXX, tft.width(),

0);

p.y = map(p.y, TS_MINY, TS_MAXY, tft.height(),

0);

int xt = p.x/cellSize;

int yt = p.y/cellSize;

vt[xt][yt] = 3;

v[xt][yt] = 3;

for(int i = 0; i < PENRADIUS; i++){

for(int j = 0; j < 20; j++){

xt = (i * cos(j*PI/10) + p.x)/cellSize;

yt = (j * sin(j*PI/10) + p.y)/cellSize;

vt[xt][yt] = 3 ;

v[xt][yt] = 3;

}

}

}

Here, the interacted cells are set to a constant value $v=3$⁠, which is comparable to the value of excited cells in the simulation. These interactions begin as target waves propagating outward. Collision with other waves can cause the spiral waves to break up, leading to more complex dynamics.

The Ardiuno simulations can also be directly compared with, for example, the spiral waves formed with the Belousov-Zhabotinsky (BZ) chemical reaction.⁶⁶ In Fig. 6, the FHN spirals are compared with those formed by the BZ chemical reaction in a Petri dish (see the supplementary material for information on how to make the BZ chemical reaction shown here). These spiral waves are produced when concentric waves activated by a silver wire, which break due to interactions with gradients in refractory periods.^67–69 Simulations of the activations by the silver wire are implemented in the Arduino by touching the screen.

The complete sketch for this example is listed in the supplementary material as “FHN_2D.ino.”

The second model implemented is the Fenton-Karma (FK) model.⁸ This model reproduces cardiac action potentials using three variables, $V$⁠, $v$⁠, and $w$⁠, and captures similar behavior to that produced by some more complex cardiac cell models like the Beeler-Rueter (BR) or Luo-Rudy (LR) I models, both of which have eight variables but without the heavier computational load. This model reproduces several key aspects of cardiac tissue, including the upstroke time scale of the AP, adaptation of the AP and conduction velocity to changes in the pacing rate, and a minimal diastolic interval before conduction block. These properties, which are set through tunable parameters, allow the model to match not only the dynamics of other cardiac cell models but also to reproduce experimental data.⁵⁸ 

The FK model is given by a set of three differential equations,

where

In this model, $Ifi$⁠, $Iso$⁠, and $Isi$ represent the fast inward $Na+$ current, the slow inward $Ca2+$ current, and the slow outward $K+$ current, respectively. The $V,v$⁠, and $w$ variables represent the membrane voltage, fast sodium gate, and slow calcium gate variables, respectively. All other parameters are described in the original publication⁸ and the programs in the supplementary material. The voltage in this model, as that in the FHN, is normalized to have amplitudes close to 1. To compare with realistic transmembrane voltage values found in cardiac cells (given in millivolts), the $V$ value can be rescaled by $V∗100mV−85mV$⁠.

One key complex dynamical instability that appears in cardiac tissue in most animal species is a period-doubling bifurcation known as alternans that develops at fast pacing periods of stimulation. This bifurcation was first described via a cobweb map⁴⁰ of the action potential restitution curve by Nolasco and Dahlen and then using linear stability analysis by Guevara et al.⁴¹ This bifurcation is important as it has been linked clinically to T-wave alternans (TWA),⁴² a prognostic for arrhythmias. People diagnosed with TWA have been shown to have a very short survival life span (80% death within 2 years).⁷⁰ TWA is also used by the U.S. Food and Drug Administration as a marker for drugs that can lead to arrhythmias. Thus, over the years, there have been many efforts toward the study of this bifurcation in time and space.^{39,42,70–72}

One example of alternans can be observed in rings when pulses change wavelength due to this bifurcation.³⁹ When the ring is small enough for the period of rotation to fall below the critical period where the period-doubling bifurcation occurs, discordant alternans develops.^23,59 Alternans is shown in Fig. 7 for the LEDs and a Java applet where successive pulses alternate between long and short periods as they propagate,⁷³ in contrast with the constant wavelength shown by the FHN. The LED ring representation gives a clear visualization of this shrinking and lengthening of pulses. Panels (a)–(c) $t1−t4$ show one pulse that grows in time as it propagates until it reaches the back of the same wave. Furthermore, at very high levels of alternans, sometimes wave backs collapse⁷³ resulting in a single wave front and three wave backs [Figs. 7(a)–7(c)] $t5−t6$⁠. This collapsing wave phenomenon has recently been shown in canine experiments.⁷⁴ For the size of the ring here, the wave ultimately collapses fully.

In the FK model, alternans can be observed for particular sets of parameters⁷⁵ as a function of the size of the ring (which can be accomplished here by changing, for example, the diffusion coefficient or the spatial resolution while keeping the integration resolved),²³ as shown in Fig. 7.

The complete sketch for this example is listed in the supplementary material as “FK.ino.”

The Kuramoto model⁷⁶ (KM) is a widely explored model used to describe the dynamics of oscillators with weak coupling that synchronize when the coupling strength is larger than a critical value.⁷⁷ It was originally motivated to describe biological and chemical oscillators,^76,78 but it can also describe the behavior of physical oscillator systems like coupled pendula and an array of coupled Josephson junctions.⁷⁹ This system can be studied with various topologies or complex networks^80–82 as well as with noise or time delay,⁸³ depending on the system of interest. Demonstrated here is the Arduino implementation of the KM for three cases: the completely connected (case 1), the spatially extended two-dimensional topology with the nearest-neighbor coupling (case 2), and the ring topology with nonlocal coupling (case 3), which exhibits a complex dynamics called “chimeras.”

Each Kuramoto oscillator obeys the following differential equation for its phase:

This model describes the change of phase $θi$ of oscillators with natural frequency $ωi$ and coupling strength $σ$⁠. The regime of interest is characterized by coupling that is small compared to the natural frequency. For this simulation, $ωi$ is sampled from the Gaussian distribution centered at $3π$ in normalized units with a standard deviation of $π/2$ and $σ$ in the interval $[0,1]$⁠. Each oscillator is coupled to every other oscillator in order to exhibit the effects of global coupling.

The critical coupling, $σc$⁠, is the coupling after which the steady-state dynamics change. For coupling larger than the critical coupling (⁠$σ>σc>0$⁠), the system will always phase-synchronize in the completely connected model. Below this critical coupling, $σc>σ>0$⁠, the oscillators will remain decoherent.

In this system, phase coherence can be studied through an order parameter^77,84 that can be defined as follows:

In this system, the oscillators can be visualized as traveling over the interval $[0,2π)$ with an average phase of $ϕ$⁠, with $r$ describing the coherence of the oscillators. When $r=0$⁠, the phases are completely spread out, and when $r=1$⁠, the phases are all the same.

For this setup, an $8×8$ LED grid and the 24-LED ring are used. In the $8×8$ grid, each LED represents a single Kuramoto oscillator inspired by thinking about the oscillator on a unit circle and mapping it to an RGB color wheel. The colors of each LED in the grid are determined as follows:

for(int i = 0; i < L; i++){

for(int j = 0; j < L; j++){

  int red = 0;

  int blue = 0;

  int green = 0;

  if(theta[i][j] <= 2 * pi/3){

   green = theta[i][j] * 255/(2 * pi);

   red = 255 - theta[i][j] * 255/(2

   * pi/3);

   blue = 0;

  }else if(theta[i][j] <= 4 * pi/3){

   green = 255 - (theta[i][j]

   - 2 * pi/3) * 255/(2 * pi/3);

   blue = (theta[i][j]

   - 2 * pi/3) * 255/(2 * pi/3);

   red = 0;

  }

  else{

   blue = 255 - (theta[i][j]

   - 4 * pi/3) * 255/(2 * pi/3);

   red = (theta[i][j]

   - 4 * pi/3) * 255/(2 * pi/3);

   green = 0;

  }

  matrix.setPixelColor(i + j * L, red,

  green, blue);

 }

}

For the 24-LED ring, the top LED corresponds to the range $[0,π12)$ moving clockwise and the intensity of the LED light is proportional to the number of oscillators in that particular phase range. To code this scenario, the following lines are included:

double ang = (theta[i][j] * 180/pi)/15;

int phase = (int) ang;

for (int j = 0; j < 24; j++){

  if (phase == j){

   phi[j]++;

}

}

In this way, the circle is divided into 24 sections, one for each LED in the ring, which is shown in Fig. 8. When the oscillator is determined to have a phase in section $i$⁠, the count $phi[i]$ is increased by one. This process bins the oscillators to be represented by a particular LED when its phase is in the range assigned to that LED. The more oscillators are in a given bin, the brighter the LED. In this way, it will be clear when the oscillators become phase-synchronized. The colors of the LEDs are determined as follows:

for (int i = 0; i < 24; i++){

  int red = 0;

  int blue = 0;

  int green = 0;

  if(i < 8){

   red = (255 - (i * 32 - 1)) * phi[i]/64;

   green = (i * 32 - 1) * phi[i]/64;

   blue = 0;

  }else if(i < 16){

   green = (255 - ((i - 8) * 32 - 1))

   * phi[i]/64;

   blue = ((i - 8) * 32 - 1) * phi[i]/64;

   red = 0;

  }else {

   blue = phi[i] * (255 - (i - 16) * 32

   - 1)/64;

   red = ((i - 16) * 32 - 1) * phi[i]/64;

   green = 0;

  }

  strip.setPixelColor(i, red, green, blue);

  if(phi[i]==0){

   strip.setPixelColor(i, 0, 0, 0);

  }

}

Each bin is displayed by a different color in a gradient between red, green, and blue. The intensity of the color denotes the number of oscillators in this bin, where $0$ represents no oscillators.

When running the simulation on the Arduino and LED ring, a stopwatch is used to time how quickly a range of different coupling constants will reach phase coherence, which is displayed in the last snapshot of Fig. 8. This time to steady state behavior is denoted as $t∞$⁠. These values are plotted as $ln⁡(σ)$ vs $ln⁡(t∞)$ in the black circles in Fig. 10 and show a linear log-log relationship.

Each Kuramoto oscillator obeys the following differential equation for its phase:

For this case, we imagine the network as a two-dimensional grid of oscillators with phases $θi,j$⁠, where $i,j=0,1,…,7$⁠, $i$ describes the row and $j$ describes the column. Periodic boundary conditions are enforced by identifying boundaries such that $i,j=−1=7$ and also $i,j=8=0$⁠.

The oscillators have natural frequencies $ωi,j$ and coupling strengths $σi,j$ with the same ranges as in case 1. In this system, the order parameter is expressed to take into account the position,

The difference between case 1 and case 2 lies in the value of the critical coupling constant $σc1$ and $σc2$⁠, which is smaller for case 1 than for case 2.

With the Adafruit $8×8$ LED matrix and the 24-LED ring, synchronization of the oscillators can be observed. Figure 9 shows Java simulations and the LEDs that are controlled by the Arduino for the nearest-neighbor coupling case. The system in this simulation takes a time equivalent to many periods of the oscillators to reach synchronization. During these intermediate steps, sections of oscillators will begin phase-synchronizing with their neighbors, which can be seen in the patches of LEDs with the same color.

The times to steady-state behavior, in this case phase coherence, $t∞$⁠, in both the completely connected and the nearest-neighbor coupling case are shown in Fig. 10. The completely connected topology reaches coherence faster than the nearest-neighbor coupling case. Also, $ln⁡(σ)∝−ln⁡(t∞)$⁠. The time to synchronization in the nearest-neighbor coupling case is longer for each value of coupling $σ$ because each oscillator has fewer connections (4 vs 63 in the completely connected case), meaning that smaller sections become synchronized en route to complete phase synchronization. Figure 9 shows that some spread in the oscillator phase remains. In some cases, splay states⁸⁵ may even emerge, where the spread is larger and ultimately is spread out more uniformly in the system, causing a periodic wave to travel through the oscillators.

The complete sketch for the previous two examples is listed in the supplementary material as “Kuramoto.ino.”

“Chimeras” are states of systems that are all composed of identical oscillators with the same dynamics, but, when coupled in space, they can develop regions with different dynamics.⁸⁶ Typically, systems that exhibit chimera states use nonglobal coupling.^86,87

For these simulations, the initial conditions and coupling parameters were taken from Abrams and Strogatz,⁸⁶ which also includes a feedback parameter $α$⁠,

The coupling is given by

where $x$ is the oscillator number and $0≤A≤1$⁠. This coupling form makes oscillators that are close more strongly coupled than those farther from each other. Initial conditions are found by a random distribution of

where $r$ is a uniform variable on $[−12,12]$⁠. This distribution yields oscillators with small fluctuations at the oscillators near the boundaries, making them closer to synchronizing, while the other oscillators are random.

Here, even though the oscillators are identical, the dynamics of the oscillators are not. In particular, there are two overall dynamics occurring simultaneously: coherence and decoherence, i.e., a type of chimera state. The 60-LED ring is used to show these dynamics in Fig. 11. In this image, the 60-LED ring represents the physical space of the oscillators and the color of the LED represents the phase. On the left side of the ring in Fig. 11 where all the LEDs are in the same color, the oscillators are coherent and the time series shows these oscillators are synchronized and oscillating at a fixed frequency. This frequency can be measured with a stopwatch which will give the real time of the system or by printing out the values in the serial port which would give the computer time. The right side of the LEDs shows oscillators that are different colors and so are decoherent and remain so while oscillating at a certain frequency. The images in Fig. 11 are taken at every 0.34 s in real time yielding a period of 1.7 s (5 images before it repeats). The phase-coherent oscillators will remain phase-coherent, and the decoherent oscillators will remain decoherent for a long integration time.

The complete sketch for this example is listed in the supplementary material as “KM_Chimeras.ino.”

In Secs. III–V, we show how MC can interact with simple external devices. However, MCs are becoming extensively used in experimental setups where they can be used as affordable acquisition boards sending/receiving digital or analog signals to various connected devices and for closed-loop control.^88–92 Similarly, HTML pages now provide a high level of interactivity through JavaScript and can run on virtually any computer, tablet, and cell phone. They can also utilize all the computing resources that are shared by the operating system with the browser. The HTML pages can even provide utilization of graphics processing units (GPUs) through WebGL.^93,94 This technology can be a powerful tool to provide a high-performance-computing platform on personal devices that turns them into virtual supercomputers. WebGL, in conjunction with MCs, has the potential to enable us to create complicated interactive control algorithms in real time. We present the first step to implement such setups, which is a communication bridge between the MC board and the HTML page, and an example using the FHN model.

HTML pages cannot communicate directly with connected devices on a machine that runs them due to security concerns. However, they can easily communicate with web servers, specifically the server that hosts the pages. The server acts as an application that can allow communication with various external connected devices, for example, to send the HTML file/application and all its necessary files to the browser of the connected machine. If the server is allowed to communicate with the hardware, it can act as a message relay demon between the MC board and the HTML page (see Fig. 12).

In this example, we will be using Node.js for the web server, as it provides bidirectional communication between the server and the HTML page through a web-socket protocol as well as between the server and the Arduino MC board through the serial port. Our objective is to provide an HTML front end for a single-cell FHN solver on the Arduino. The front end will be able to visualize the voltage signal that is calculated on the MC in the browser in real time and set a range of modifiable parameters on the MC. Therefore, the Node.js server should be able to handle at least three tasks: (1) serving the HTML file and its dependencies to the browser, (2) receiving messages from the HTML page through socket.io and relaying them to the Arduino board through the serial port, and (3) receiving messages through the serial port from the MC and relaying this information to the HTML program to display. The following code is the Node.js setup file used to achieve all these steps; it includes ample comments for clarity.

var express = require(’express’);

var app = express() ;

var http = require(’http’).Server(app);

app.get(’/’, function(req, res){

res.sendFile(__dirname+’/www/index.html’);

});

app.use(express.static(’www’)) ;

http.listen(3000, function(){

console.log(’listening on *:3000’);

});

/* Importing serialport package as

SerialPort */

var SerialPort = require(’serialport’) ;

const Readline = require(’@serialport/parser

- readline’) ;

const parser = new Readline() ; /* parses

messages from the Arduino */

/* Serial port address that the MC board is

connected to. It can be with the format

’COM3’ on Microsoft Windows, or

’/dev/ttUSB1’ on Linux, or

’/dev/cu.usbserial-1420’, on Mac OS. */

var port_address = ’/dev/cu.usbserial-1420’ ;

/* Defining the serial port for bi-directional

communication between the server and the MC.*/

var port = new SerialPort(port_address,{

 baudRate: 9600,

 dataBits: 8,

 parity : ’none’,

 stopBits: 1,

 flowControl : false,

}) ;

port.pipe(parser) ;

/* The parser will send each line of data

that is received through the serial port

from the MC to the HTML page through the

socket.io with the flag ’A’. */

parser.on( ’data’, line =>

io.emit(’A’,‘{line}‘)) ;

/* Setting up the socket.io to set

bi-directional communication between the

server and the HTML application. */

var io = require(’socket.io’)(http) ;

  io.on(’connection’, function(socket){

  console.log(’a user connection is

  established’) ;

  socket.on(’disconnect’, function(){

   console.log(’a user

   disconnected’) ;

  } ) ;

  socket.on(’msg’, function(msg){

   /* When a message is received from

   the HTML application, relay the

   message to the MC through the

   serial port */

   port.write(msg) ;

   });

} ) ;

The HTML application will require the inclusion of the socket.io JavaScript library. The following line of code imports that library.

<script src="/socket.io/socket.io.js"></script>.

Then, a bidirectional communication channel is defined through JavaScript in the HTML file achieved by the following line:

var socket = io() ;

In this example, the only information that is received from the MC board is the membrane potential, which is flagged by the string “A.” Therefore, each time a voltage value is received, time is incremented by the discretization time step and the voltage plot is updated with the received value of voltage that is achieved by the following JavaScript code in the HTML application. Additionally, the time access is reset every “finalTime” units of time to be able to continuously plot the voltage signal.

/* This snippet defines a message-received

event:

When the message (msg) with the flag ’A’

is received, the function(msg){...} is

executed. */

socket.on(’A’, function(msg){

  time += dt ; /* total time */

  ptime += dt ; /* plot time */

  /* The variable msg is a String and

  needs to be parsed as a float to

  be plotted. */

  ucrv.plot(time,parseFloat(msg)) ;

  /* Each time, ptime passes the

   finalTime, the plotting window is

   slided by finalTime to have a

   continuous plot of the voltage. */

   if(ptime > finalTime){

   ptime = 0 ;

   plt.xlimits = [ plt.xlimits[0]

   + finalTime,plt.xlimits[1]

   + finalTime ] ;

}

} ) ;

To create an interactive graphical interface, regular input elements such as buttons and numbers can be added to the page. Whenever an interaction with the page requires informing the MC about the interaction, a message with the appropriate code is sent to the server to be relayed to the MC. For example, in the following snippet, when an element with the id “a” changes, a message with the flag “ms” and starting with “a” followed by the value of the element is sent to the server to be relayed to the MC.

document.getElementById(’a’).onchange

= function(){

 socket.emit( ’msg’, ’a’+this.value+’&’ ) ;

} ;

It should be noted that the MC reads the serial data one character at a time; hence, even floating-point numbers are sent as a series of characters. Therefore, when we are done writing a number, we add a “&” character to mark the end of the number or message to facilitate message parsing on the MC. To achieve this, the following function is implemented.

float readnumber(){

 String incoming ;

 char ch ;

 while ((ch = Serial.read())!=’&’ ){

  if( isDigit(ch) || ch == ’.’ ){

   incoming += (char) ch ;

  }

 }

  return incoming.toFloat() ;

}

A switch statement is used to parse the data that is sent to the MC from the server. A simple graphical interface is designed to interact with the MC board through the web page (see Fig. 13). To showcase the ease of use of the setup, different pacing periods were tested with the parameters shown in Fig. 13. Figure 14 shows this example. As shown, the action potential duration (APD) becomes shorter and shorter until alternans (long-short APDs) develops⁴¹ as the pacing period is shortened in each panel of the figure.

While the example shown here is a relatively straightforward, showing how a microcontroller can interact in real time with a web-browser opens up a large range of possibilities for simulations of devices with open-loop feedback control using fast WebGL simulations,^93,94 a microcontroller, and a physical external device.

We have presented here a way to implement and study complex systems, including reaction-diffusion systems in one and two dimensions, using easy and inexpensive microcontrollers connected with LEDs, touchscreens, and a web browser. These tools open the door to alternative, more interactive ways to study and teach complex systems. It especially allows for activities with a large number of students, each with their own system or a few working on the same system as we have done in several national and international workshops or with undergraduate summer research students. With the codes presented here, students can easily investigate traveling waves, spiral waves, and coupled-oscillator synchronization. These codes can be used in a laboratory setting where students change parameters such as the coupling strength to study how they affect the synchronization for different scenarios, such as a 2D array of oscillators with the nearest-neighbor coupling in both physical and phase space, through the Adafruit matrix, ring, and LCD displays. In addition, other more complex models can easily be implemented using more advanced microcontrollers.

Students can expand on this setup and are also able to change how the oscillators are coupled in the Arduino sketches for the different models implemented as well as to mix the examples up such that the FK or KM model includes interactions via light or touch, which we did not include here. Furthermore, using other sorts of LED setups, students can create their own networks of oscillators using the same code to simulate the KM, even making unilateral coupling or coupling with different strengths. Also, the 2D KM can be written to work with the touchscreen if desired, and the screens can be set up to better display other sorts of topologies.

Some open projects not presented here include the use of an LED light to be measured by a photoresistor and with coupling to neighboring LED oscillators dependent on the value measured by the photoresistor, rather than coupling determined from a set diffusion value, which can allow external real random noise to be part of the dynamics. Another example could include the use of the MC connected to the web page to run a highly complex simulation in WebGL in real time using the Arduino as a controller for feedback. For example, the FK model can be run during alternans in WebGL to simulate waves propagating in a cardiac Purkinje fiber, with the system running in real time at 200 cm/s and producing alternans,⁹⁵ and the Arduino can be used to perform control in the loop.^96–100 

See the supplementary material for detailed explanation on Arduino connections, BZ reaction, and access to all the Arduino sketches discussed in this paper.

This work was supported, in part, by the National Science Foundation (NSF) (No. CNS-1446675 to F.H.F.) and by the National Institutes of Health (NIH) (Grant No. 1R01HL143450-01 to F.H.F.). F.H.F. would also like to thank the NIMBioS Working Group on “Prediction and control of cardiac alternans” (see Ref. 101) for useful discussions.