The spatiotemporal dynamics of complex systems have been studied traditionally and visualized numerically using high-end computers. However, due to advances in microcontrollers, it is now possible to run what once were considered large-scale simulations using a very small and inexpensive single integrated circuit that can furthermore send and receive information to and from the outside world in real time. In this paper, we show how microcontrollers can be used to perform simulations of nonlinear ordinary differential equations with spatial coupling and to visualize their dynamics using arrays of light-emitting diodes and/or touchscreens. We demonstrate these abilities using three different models: two reaction-diffusion models (one neural and one cardiac) and a generic model of network oscillators. These models are commonly used to simulate various phenomena in biophysical systems, including bifurcations, waves, chaos, and synchronization. We also demonstrate how simple it is to integrate real-time user interaction with the simulations by showing examples with a light sensor, touchscreen, and web browser.

Microcontrollers (MCs) are single integrated chips that contain a powerful central processing unit (CPU), memory, a clock, and many input-output control units. In essence, they are fast computers on single circuits that can digitally control other processes and are increasingly integrated into cyberphysical systems from automobiles to toys. Due to their speed and affordability relative to traditional computers, they are also used in scientific computation. Examples include the logistic map1 and closed-loop experiments of interacting waves in the Belousov-Zhabotinsky (BZ)2 chemical reaction. More recently, they have been used to simulate electrical waves in one-dimensional cardiac tissue in real time3,4 by connecting an array of microcontrollers in series, where each microcontroller represents a single cardiac cell.5 However, the use of microcontrollers for applications is still in its early stages and much remains to be explored, particularly regarding the use of emerging peripherals and readout interfaces. In this paper, we connect a microcontroller to an array of light-emitting diodes (LEDs) and a liquid crystal display (LCD) to simulate and visualize the dynamics of three models: (1) the FitzHugh-Nagumo (FHN) model,6,7 (2) the three-variable Fenton-Karma (FK) model of cardiac dynamics,8 and (3) a 2D spatial network of coupled Kuramoto oscillators. The FHN and FK models are visualized as waves in a ring, while the Kuramoto model (KM) synchronization phenomenon is visualized on both a spatial grid and a phase-space representation. Furthermore, we demonstrate the ease of direct interaction between the real world and the simulation by implementing sensors as well as interaction with a web browser that allows for real-time manipulation of the system’s state during dynamical evolution. We present not only a new, simple, and scalable method to study and visualize complex systems with and without the use of computers and monitors but also a means of coupling these complex simulated systems to the local environment of the controller itself.

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.15 These systems can span from simple models such as FHN6,7 and Hodgkin-Huxley (HH),16 consisting of two and four variables, respectively, to complex models with more than 80 variables.17 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,21 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.27 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.1 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.1 Separately, Mahmud et al.3 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.30 MCs also have been used to simulate the HH neuron model with output via the serial port or to an LCD;31 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.37 discuss the use of different MCs in robotics for classroom activities and in the advancement of automatic control.37 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.38 

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.16 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.8 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 alternans39 comparable to those observed in cardiac experiments.40–43 

Next, we use MCs to illustrate synchronization,44 a spontaneous emergence of order that can be observed in coupled oscillators such as some cardiac cells45 and fireflies.46 We use the Kuramoto model,47 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 seizures48 and the resonance caused by people crossing the London Millennium Footbridge while inadvertently synchronizing their footsteps with the oscillation of the structure.49 

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.50 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.

FIG. 1.

Arduino Due connected to the Adafruit NeoPixel 24-RGB LED ring and the NeoPixel 8×8 RGB LED matrix. Other setups are shown and discussed in detail in the supplementary material.

FIG. 1.

Arduino Due connected to the Adafruit NeoPixel 24-RGB LED ring and the NeoPixel 8×8 RGB LED matrix. Other setups are shown and discussed in detail in the supplementary material.

Close modal

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,51 neurons,52 and cardiomyocytes.53 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,

vt=v(av)(v1)w,wt=ϵ(bvdwδ).
(1)

Parameters a, b, d, δ, and ϵ affect the location and type of equilibria and the speed of oscillations.6 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.

FIG. 2.

Six images of the touchscreen showing the action potential of one FHN oscillator. (a) A few seconds after the simulation begins, a full action potential is drawn on the touchscreen. (b) This action potential will repeat at a constant frequency if left unperturbed (untouched). (c) While interacting with the touchscreen, the oscillator will excite (when touched), causing the voltage to increase by a certain amount depending on the pressure applied to the screen. (d) Releasing the touchscreen will cause the oscillator to start to relax. When the oscillator is excited too soon after the previous excitation, the action potential will be shorter than the typical size. (e) and (f) Fast perturbations by interacting with the touchscreen very quickly after the oscillator relaxes will cause smaller and shorter excitations as the oscillator is perturbed closer to the refractory period. The parameters used for this simulations are a=2.5, b=2.5, ϵ=1, d=0, and δ=0 with dt=0.05.

FIG. 2.

Six images of the touchscreen showing the action potential of one FHN oscillator. (a) A few seconds after the simulation begins, a full action potential is drawn on the touchscreen. (b) This action potential will repeat at a constant frequency if left unperturbed (untouched). (c) While interacting with the touchscreen, the oscillator will excite (when touched), causing the voltage to increase by a certain amount depending on the pressure applied to the screen. (d) Releasing the touchscreen will cause the oscillator to start to relax. When the oscillator is excited too soon after the previous excitation, the action potential will be shorter than the typical size. (e) and (f) Fast perturbations by interacting with the touchscreen very quickly after the oscillator relaxes will cause smaller and shorter excitations as the oscillator is perturbed closer to the refractory period. The parameters used for this simulations are a=2.5, b=2.5, ϵ=1, d=0, and δ=0 with dt=0.05.

Close modal

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:

vt=v(av)(v1)w+I,wt=ϵ(bvdwδ),
(2)

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.6 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 exist54 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.57 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,58 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.59 In the FHN model, there is a small parameter regime where waves can be elicited intermittently in each direction60 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,

vit=vi(avi)(vi1)w+D2vi,wit=ϵ(bvidwiδ).
(3)

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.61 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.

FIG. 3.

(a) Six images at different times of the Adafruit 60-LED ring running FHN model, where the brightest red LEDs show the peak v and the off LEDs show the minimum v. (b) Six still images of the simulation of the FHN model on a ring taken at 25 time step intervals. (c) Voltage-position plots show that the wave propagates clockwise and the wavelength is constant.

FIG. 3.

(a) Six images at different times of the Adafruit 60-LED ring running FHN model, where the brightest red LEDs show the peak v and the off LEDs show the minimum v. (b) Six still images of the simulation of the FHN model on a ring taken at 25 time step intervals. (c) Voltage-position plots show that the wave propagates clockwise and the wavelength is constant.

Close modal

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.61 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.62 

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)].

FIG. 4.

Interaction with a photocell, shown in the blue circle, with a FHN wave that is traveling around the LED ring. (a) The FHN wave is traveling counterclockwise. (b) Light to the photocell is blocked via a finger and cells near the photocell light up as they become activated. (c) This activation begins to spread as a counterclockwise wave, while the cells in the clockwise direction are still in the refractory period and block the wave. (d) The initial wave and the new wave collide. (e) The two waves begin to annihilate each other. (f) None of the cells are activated, and the traveling wave is gone.

FIG. 4.

Interaction with a photocell, shown in the blue circle, with a FHN wave that is traveling around the LED ring. (a) The FHN wave is traveling counterclockwise. (b) Light to the photocell is blocked via a finger and cells near the photocell light up as they become activated. (c) This activation begins to spread as a counterclockwise wave, while the cells in the clockwise direction are still in the refractory period and block the wave. (d) The initial wave and the new wave collide. (e) The two waves begin to annihilate each other. (f) None of the cells are activated, and the traveling wave is gone.

Close modal

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)].

FIG. 5.

Example 2 of interaction with a photocell, shown in the blue circle, with FHN oscillators. (a) There is no initial wave. Light to the photocell is blocked via a finger, and cells near the photocell light up as they become excited. (b) and (c) This activation begins to spread as a wave both clockwise and counterclockwise. (d) and (e) The two wave fronts collide, and the wave begins to die out. (f) The wave dies out completely, and there is no more activation.

FIG. 5.

Example 2 of interaction with a photocell, shown in the blue circle, with FHN oscillators. (a) There is no initial wave. Light to the photocell is blocked via a finger, and cells near the photocell light up as they become excited. (b) and (c) This activation begins to spread as a wave both clockwise and counterclockwise. (d) and (e) The two wave fronts collide, and the wave begins to die out. (f) The wave dies out completely, and there is no more activation.

Close modal

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.63 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 fibrillation64 and retinal spreading depression.65 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.61 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.63 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.66 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.

FIG. 6.

Top: Images of the 2D FHN model running on the Arduino with the LCD. (a) Initial conditions used to begin a spiral wave. Half a row begins as refractory, shown as the darker color, with a row above slightly excited (not seen). All other cells are at rest. (b) A user interacting with the display as the simulation continues. (c) Three white spots on the right of the display show three locations where the user interacted with the display; those cells are now excited. (d) The three locations propagate as target waves outward from the interacted cells. The initial activation beings to propagate out and rotates around the initial refractory cells. (e) Two spirals begin to form as the previous waves collide. (f) Much later after two stable spiral waves have formed. Bottom: Spiral waves in the Belousov-Zhabotinsky chemical reaction. The white color is the activation in the reaction, and the waves are able to propagate as spirals.

FIG. 6.

Top: Images of the 2D FHN model running on the Arduino with the LCD. (a) Initial conditions used to begin a spiral wave. Half a row begins as refractory, shown as the darker color, with a row above slightly excited (not seen). All other cells are at rest. (b) A user interacting with the display as the simulation continues. (c) Three white spots on the right of the display show three locations where the user interacted with the display; those cells are now excited. (d) The three locations propagate as target waves outward from the interacted cells. The initial activation beings to propagate out and rotates around the initial refractory cells. (e) Two spirals begin to form as the previous waves collide. (f) Much later after two stable spiral waves have formed. Bottom: Spiral waves in the Belousov-Zhabotinsky chemical reaction. The white color is the activation in the reaction, and the waves are able to propagate as spirals.

Close modal

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.8 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.58 

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

tV(x,t)=(D~V)Ifi(V;v)+Iso(V)+Isi(V;w)Cm,tv(x,t)=(1Θ(VVc))(1v)τv(V)Θ(VVc)vτv+,tw(x,t)=(1Θ(VVc))(1w)τw(V)Θ(VVc)wτw+,
(4)

where

Ifi(V;v)=vΘ(VVc)(VVc)(VmV)τd,Iso(V)=(VVo)1Θ(VVc)τo+Θ(VVc)τr,Isi(V;w)=w1+tanh(k(VVcsi))2τsi.
(5)

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 publication8 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 V100mV85mV.

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 map40 of the action potential restitution curve by Nolasco and Dahlen and then using linear stability analysis by Guevara et al.41 This bifurcation is important as it has been linked clinically to T-wave alternans (TWA),42 a prognostic for arrhythmias. People diagnosed with TWA have been shown to have a very short survival life span (80% death within 2 years).70 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.39 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,73 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) t1t4 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 collapse73 resulting in a single wave front and three wave backs [Figs. 7(a)7(c)] t5t6. This collapsing wave phenomenon has recently been shown in canine experiments.74 For the size of the ring here, the wave ultimately collapses fully.

FIG. 7.

(a) A column of six still images taken at 48-s intervals of the 60-LED ring simulating the FK model. The red LEDs are the cells that are activated. The wavefront travels counterclockwise. (b) A column of six frames from a computer simulation of the FK model on a ring with images taken at intervals of 48 time steps. The white sections of the ring are the parts that are activated. (c) The voltage-position plots of the FK simulations. The wave propagates counterclockwise and the wavelength oscillates in time. The 5th set of images shows “alternans” as the development of wavelength oscillations; here, the wave back collapses in the middle and forms two additional wave backs. The 6th set shows one wave dying out and the other continuing to travel counterclockwise.

FIG. 7.

(a) A column of six still images taken at 48-s intervals of the 60-LED ring simulating the FK model. The red LEDs are the cells that are activated. The wavefront travels counterclockwise. (b) A column of six frames from a computer simulation of the FK model on a ring with images taken at intervals of 48 time steps. The white sections of the ring are the parts that are activated. (c) The voltage-position plots of the FK simulations. The wave propagates counterclockwise and the wavelength oscillates in time. The 5th set of images shows “alternans” as the development of wavelength oscillations; here, the wave back collapses in the middle and forms two additional wave backs. The 6th set shows one wave dying out and the other continuing to travel counterclockwise.

Close modal

In the FK model, alternans can be observed for particular sets of parameters75 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),23 as shown in Fig. 7.

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

The Kuramoto model76 (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.77 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.79 This system can be studied with various topologies or complex networks80–82 as well as with noise or time delay,83 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:

dθidt=ωi+σj=0N1sin(θjθi).
(6)

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 parameter77,84 that can be defined as follows:

reiϕ=1Nj=0n1eiθj.
(7)

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:

FIG. 8.

Top: Java applet of the Kuramoto model on an 8×8 LED grid with order parameter shown by the ring. The individual oscillators traverse 0 to 2π ring, from red to green to blue-back to red. Bottom: Same setup on the Arduino 8×8 LED array and 24-LED ring. Both use global coupling. The Arduino images are taken every 4 frames from a 59 frame/s video.

FIG. 8.

Top: Java applet of the Kuramoto model on an 8×8 LED grid with order parameter shown by the ring. The individual oscillators traverse 0 to 2π ring, from red to green to blue-back to red. Bottom: Same setup on the Arduino 8×8 LED array and 24-LED ring. Both use global coupling. The Arduino images are taken every 4 frames from a 59 frame/s video.

Close modal

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:

dθi,jdt=ωi,j+σ[sin(θi,j+1θi,j)+sin(θi,j1θi,j)+sin(θi+1,jθi,j)+sin(θi1,jθi,j)].
(8)

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,

reiϕ=1Nj,k=0n1eiθj,k.
(9)

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.

FIG. 9.

Top: Java applet of the Kuramoto model on an 8×8 LED grid with order parameter shown by the ring. The individual oscillators traverse 0 to 2π, from red to green to blue-back to red. Bottom: Same setup on the Arduino 8×8 LED array and 24-LED ring. The simulations are run with nearest-neighbor coupling.

FIG. 9.

Top: Java applet of the Kuramoto model on an 8×8 LED grid with order parameter shown by the ring. The individual oscillators traverse 0 to 2π, from red to green to blue-back to red. Bottom: Same setup on the Arduino 8×8 LED array and 24-LED ring. The simulations are run with nearest-neighbor coupling.

Close modal
FIG. 10.

Data of time to synchronization of the LEDs for different coupling strengths for two coupling topologies. There is a shift in time to synchronization from the completely connected case to the nearest-neighbors case.

FIG. 10.

Data of time to synchronization of the LEDs for different coupling strengths for two coupling topologies. There is a shift in time to synchronization from the completely connected case to the nearest-neighbors case.

Close modal

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 states85 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.86 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,86 which also includes a feedback parameter α,

dθdt=ωππG(xx)sin[(θ(x,t)θ(x,t)+α]dx.
(10)

The coupling is given by

G(x)=12π(1+Acosx),
(11)

where x is the oscillator number and 0A1. 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

θ(x)=6re0.76x2,
(12)

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.

FIG. 11.

A ring of 60-LEDs running a simulation for chimera dynamics. The color of each LED represents a phase. The LEDs that are in the dashed semiarcs represent Kuramoto oscillators that remain decoherent for long periods of time but that are still oscillating in time. The LEDs outside of that arc are Kuramoto oscillators that are phase-synchronized. Panels (a)–(f) show frames spaced every 0.34 s of real time showing one complete period of 1.7 s.

FIG. 11.

A ring of 60-LEDs running a simulation for chimera dynamics. The color of each LED represents a phase. The LEDs that are in the dashed semiarcs represent Kuramoto oscillators that remain decoherent for long periods of time but that are still oscillating in time. The LEDs outside of that arc are Kuramoto oscillators that are phase-synchronized. Panels (a)–(f) show frames spaced every 0.34 s of real time showing one complete period of 1.7 s.

Close modal

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

In Secs. IIIV, 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).

FIG. 12.

The Arduino MC and the HTML application can communicate only through the Node.js server.

FIG. 12.

The Arduino MC and the HTML application can communicate only through the Node.js server.

Close modal

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) develops41 as the pacing period is shortened in each panel of the figure.

FIG. 13.

Graphical interface of the webpage interface that was implemented to interact with the MC board. The buttons can initialize the solution, start/stop the simulation, and stimulate the cell. Through the interface, users can set the parameter values in the FHN model or choose a constant period stimulation protocol.

FIG. 13.

Graphical interface of the webpage interface that was implemented to interact with the MC board. The buttons can initialize the solution, start/stop the simulation, and stimulate the cell. Through the interface, users can set the parameter values in the FHN model or choose a constant period stimulation protocol.

Close modal
FIG. 14.

Four constant pacing periods were tested with the FHN single-cell model: (a) 300, (b) 150, (c) 100, and (d) 75. The signal represents the variable u in the FHN model. Panels (a)–(c) show a 1:1 response, while panel (d) shows a 2:2 period-doubling response.

FIG. 14.

Four constant pacing periods were tested with the FHN single-cell model: (a) 300, (b) 150, (c) 100, and (d) 75. The signal represents the variable u in the FHN model. Panels (a)–(c) show a 1:1 response, while panel (d) shows a 2:2 period-doubling response.

Close modal

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,95 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.

1.
J. D.
Serna
and
A.
Joshi
,
Phys. Educ.
47
,
736
(
2012
).
2.
H.
Yokoi
,
A.
Adamatzky
,
B.
De Lacy Costello
, and
C.
Melhuish
,
Int. J. Bifurcat. Chaos
14
,
3347
(
2004
).
3.
F.
Mahmud
,
N.
Shiozawa
,
M.
Makikawa
, and
T.
Nomura
,
Chaos
21
,
023121
(
2011
).
4.
F.
Mahmud
, in 2012 IEEE EMBS Conference on Biomedical Engineering and Sciences (IECBES) (IEEE, 2012), pp. 321–325.
5.
F. H.
Fenton
and
E. M.
Cherry
,
Scholarpedia
3
,
1868
(
2008
).
6.
R.
FitzHugh
,
Biophys. J.
1
,
445
466
(
1961
).
7.
J.
Nagumo
,
S.
Arimoto
, and
S.
Yoshizawa
,
Proc. IRE
50
,
2061
(
1962
).
8.
F.
Fenton
and
A.
Karma
,
Chaos
8
,
20
(
1998
).
9.
A.
Koch
and
H.
Meinhardt
,
Rev. Mod. Phys.
66
,
1481
(
1994
).
10.
S.
Kondo
and
T.
Miura
,
Science
329
,
1616
(
2010
).
11.
A.
Deutsch
and
S.
Dormann
,
Mathematical Modeling of Biological Pattern Formation
(
Springer
,
2005
).
12.
V. S.
Zykov
and
A. T.
Winfree
,
Simulation of Wave Processes in Excitable Media
(
John Wiley & Sons, Inc.
,
1992
).
13.
J. M.
Greenberg
and
S.
Hastings
,
SIAM J. Appl. Math.
34
,
515
(
1978
).
14.
A. V.
Holden
,
M.
Markus
, and
H. G.
Othmer
,
Nonlinear Wave Processes in Excitable Media
(
Springer
,
2013
), Vol. 244.
15.
Y.
Sekerci
and
S.
Petrovskii
,
Computation
6
,
59
(
2018
).
16.
A. L.
Hodgkin
and
A. F.
Huxley
,
J. Physiol.
117
,
500
(
1952
).
17.
S. N.
Flaim
,
W. R.
Giles
, and
A. D.
McCulloch
,
Am. J. Physiol. Heart Circ. Physiol.
291
,
H2617
(
2006
).
18.
M.
Courtemanche
and
A. T.
Winfree
,
Int. J. Bifurcat. Chaos
1
,
431
(
1991
).
19.
D. F.
Richards
,
J. N.
Glosli
,
E. W.
Draeger
,
A. A.
Mirin
,
B.
Chan
,
J.-l.
Fattebert
,
W. D.
Krauss
,
T.
Oppelstrup
,
C. J.
Butler
,
J. A.
Gunnels
et al.,
Comput. Methods Biomech. Biomed. Eng.
16
,
802
(
2013
).
20.
M.
Hoffman
,
N.
LaVigne
,
S.
Scorse
,
F.
Fenton
, and
E.
Cherry
,
Chaos
26
,
013107
(
2016
).
21.
G. E.
Moore
,
Proc. IEEE
86
,
82
(
1998
).
22.
M.
Dowle
,
R.
Martin Mantel
, and
D.
Barkley
,
Int. J. Bifurcat. Chaos
7
,
2529
(
1997
).
23.
F. H.
Fenton
,
E. M.
Cherry
,
H. M.
Hastings
, and
S. J.
Evans
,
Biosystems
64
,
73
(
2002
).
24.
E.
Bartocci
,
R.
Singh
,
F. B.
von Stein
,
A.
Amedome
,
A. J. J.
Caceres
,
J.
Castillo
,
E.
Closser
,
G.
Deards
,
A.
Goltsev
,
R. S.
Ines
et al.,
Adv. Physiol. Educ.
35
,
427
(
2011
).
25.
A.
Kaboudian
,
E.
Cherry
, and
F.
Fenton
,
Sci. Adv.
5
(
3
),
eaav6019
(
2019
).
26.
R. S.
Oliveira
,
B. M.
Rocha
,
R. M.
Amorim
,
F. O.
Campos
,
W.
Meira
,
E. M.
Toledo
, and
R. W.
dos Santos
, in International Conference on Parallel Processing and Applied Mathematics (Springer, 2011), pp. 111–120.
27.
S.
Augarten
, ”The most widely used computer on a chip: The TMS 1000,” in State of the Art (
Ticknor & Fields
,
1981
).
28.
D. K.
Fisher
and
P. J.
Gould
,
Modern Instrum.
1
,
8
(
2012
).
29.
E.
Koutroulis
,
K.
Kalaitzakis
, and
N. C.
Voulgaris
,
IEEE Trans. Power Electron.
16
,
46
(
2001
).
30.
C.-H.
Luo
and
Y.
Rudy
,
Circ. Res.
68
,
1501
(
1991
).
31.
Y.
Isler
,
M.
Kuntalp
, and
G.
Gonel
, in 2009 14th National Biomedical Engineering Meeting (
IEEE
,
2009
), pp. 1–4.
32.
S.
Yaghini Bonabi
,
H.
Asgharian
,
S.
Safari
, and
M.
Nili Ahmadabadi
,
Front. Neurosci.
8
,
379
(
2014
).
33.
M.
Lu
,
J.-L.
Wang
,
J.
Wen
, and
X.-W.
Dong
, in 2016 Asia-Pacific International Symposium on Electromagnetic Compatibility (APEMC) (IEEE, 2016), Vol. 1, pp. 1115–1117.
34.
Š.
Kubínová
and
J.
Šlégr
,
Phys. Educ.
50
,
472
(
2015
).
35.
C.
Galeriu
,
S.
Edwards
, and
G.
Esper
,
Phys. Teacher
52
,
157
(
2014
).
36.
W.-H.
Kuan
,
C.-H.
Tseng
,
S.
Chen
, and
C.-C.
Wong
,
J. Sci. Educ. Technol.
25
,
427
(
2016
).
37.
A.
Soriano
,
L.
Marín
,
M.
Vallés
,
A.
Valera
, and
P.
Albertos
,
IFAC Proc. Vol.
47
,
9044
(
2014
), 19th IFAC World Congress.
38.
B.
Land
, “ECE5760: Advanced microcontroller design and system-on-chip,” 2019.
39.
M. A.
Watanabe
,
F. H.
Fenton
,
S. J.
Evans
,
H. M.
Hastings
, and
A.
Karma
,
J. Cardiovasc. Electrophysiol.
12
,
196
(
2001
).
40.
J.
Nolasco
and
R. W.
Dahlen
,
J. Appl. Physiol.
25
,
191
(
1968
).
41.
M.
Guevara
,
G.
Ward
,
A.
Shrier
, and
L.
Glass
,
IEEE Comput. Cardiol.
562
,
167
(
1984
).
42.
J. M.
Pastore
,
S. D.
Girouard
,
K. R.
Laurita
,
F. G.
Akar
, and
D. S.
Rosenbaum
,
Circulation
99
,
1385
(
1999
).
43.
I.
Uzelac
,
Y. C.
Ji
,
D.
Hornung
,
J.
Schröder-Scheteling
,
S.
Luther
,
R. A.
Gray
,
E. M.
Cherry
, and
F. H.
Fenton
,
Front. Physiol.
8
,
819
(
2017
).
44.
A.
Pikovsky
,
M.
Rosenblum
,
J.
Kurths
, and
J.
Kurths
,
Synchronization: A Universal Concept in Nonlinear Sciences
(
Cambridge University Press
,
2003
), Vol. 12.
45.
Y. C.
Ji
,
I.
Uzelac
,
N.
Otani
,
S.
Luther
,
R. F.
Gilmour
, Jr.,
E. M.
Cherry
, and
F. H.
Fenton
,
Heart Rhythm
14
,
1254
(
2017
).
46.
R. E.
Mirollo
and
S. H.
Strogatz
,
SIAM J. Appl. Math.
50
,
1645
(
1990
).
47.
International Symposium on Mathematical Problems in Theoretical Physics, Lecture Notes in Physics, edited by H. Araki, J. Ehlers, K. Hepp, H. A. Weidenmüller, and W. Beiglböck (Springer, Berlin, 1975), Vol. 39.
48.
L. G.
Dominguez
,
R. A.
Wennberg
,
W.
Gaetz
,
D.
Cheyne
,
O. C.
Snead
, and
J. L. P.
Velazquez
,
J. Neurosci.
25
,
8077
(
2005
).
49.
K.
Sanderson
, “
Millennium bridge wobble explained
,” Nat. News (2008).
50.
S.
Heath
,
Embedded Systems Design
(
Newnes
,
2003
), pp. 11–12.
51.
A. L.
Lin
,
A.
Hagberg
,
E.
Meron
, and
H. L.
Swinney
,
Phys. Rev. E
69
,
066217
(
2004
).
52.
K.
Aihara
,
T.
Takabe
, and
M.
Toyoda
,
Phys. Lett. A
144
,
333
(
1990
).
53.
B. Y.
Kogan
,
W. J.
Karplus
,
B. S.
Billett
,
A. T.
Pang
,
H. S.
Karagueuzian
, and
S. S.
Khan
,
Physica D
50
,
327
(
1991
).
54.
S. P.
Hastings
,
Q. J. Math.
27
,
123
(
1976
).
55.
G.
Klaasen
and
W.
Troy
,
SIAM J. Appl. Math.
41
,
145
(
1981
).
56.
C. K. R. T.
Jones
,
Trans. Am. Math. Soc.
286
,
431
(
1984
).
57.
J. G.
Alford
and
G.
Auchmuty
,
J. Math. Biol.
53
,
797
(
2006
).
58.
E. M.
Cherry
and
F. H.
Fenton
,
New J. Phys.
10
,
125016
(
2008
).
59.
60.
E. N.
Cytrynbaum
and
T. J.
Lewis
,
SIAM J. Appl. Dyn. Syst.
8
,
348
(
2009
).
61.
Y. C.
Ji
and
F. H.
Fenton
,
Am. J. Phys.
84
,
626
(
2016
).
62.
ladyada, “Photocells,” 2012, see https://learn.adafruit.com/photocells/overview; accessed 15 January 2019.
63.
A. J.
Welsh
,
E. F.
Greco
, and
F. H.
Fenton
,
Phys. Today
70
,
78
(
2017
).
64.
J. M.
Davidenko
,
A. V.
Pertsov
,
R.
Salomonsz
,
W.
Baxter
, and
J.
Jalife
,
Nature
355
,
349
(
1992
).
65.
N.
Gorelova
and
J.
Bureš
,
J. Neurobiol.
14
,
353
(
1983
).
66.
A. M.
Zhabotinsky
,
Biofizika
9
,
306
(
1964
).
67.
A. M.
Zhabotinsky
and
A. N.
Zaikin
,
Oscillatory Processes in Biological and Chemical Systems II
(Science Publ., Puschino, 1971).
68.
A. M.
Zhabotinsky
and
A. N.
Zaikin
,
J. Theor. Biol.
40
,
45
(
1973
).
69.
A. T.
Winfree
,
Science
175
,
634
(
1972
).
70.
D. S.
Rosenbaum
,
L. E.
Jackson
,
J. M.
Smith
,
H.
Garan
,
J. N.
Ruskin
, and
R. J.
Cohen
,
New Engl. J. Med.
330
,
235
(
1994
).
71.
A.
Gizzi
,
E.
Cherry
,
R. F.
Gilmour
, Jr.,
S.
Luther
,
S.
Filippi
, and
F. H.
Fenton
,
Front. Physiol.
4
,
71
(
2013
).
72.
M.
Minkkinen
,
M.
Kähönen
,
J.
Viik
,
K.
Nikus
,
T.
Lehtimäki
,
R.
Lehtinen
,
T.
Kööbi
,
V.
Turjanmaa
,
W.
Kaiser
,
R. L.
Verrier
et al.,
J. Cardiovasc. Electrophysiol.
20
,
408
(
2009
).
73.
F. H.
Fenton
,
E. M.
Cherry
, and
L.
Glass
,
Scholarpedia
3
,
1665
(
2008
).
74.
S.
Filippi
,
A.
Gizzi
,
C.
Cherubini
,
S.
Luther
, and
F. H.
Fenton
,
Europace
16
,
424
(
2014
).
75.
F. H.
Fenton
,
E. M.
Cherry
,
H. M.
Hastings
, and
S. J.
Evans
,
Chaos
12
,
852
(
2002
).
76.
Y.
Kuramoto
, in International Symposium on Mathematical Problems in Theoretical Physics (Springer, 1975), pp. 420–422.
77.
J. A.
Acebrón
,
L. L.
Bonilla
,
C. J.
Pérez Vicente
,
F.
Ritort
, and
R.
Spigler
,
Rev. Mod. Phys.
77
,
137
(
2005
).
78.
Y.
Kuromoto
,
Chemical Oscillations, Waves, and Turbulence
(
Springer
,
1984
).
79.
K.
Wiesenfeld
,
P.
Colet
, and
S. H.
Strogatz
,
Phys. Rev. E
57
,
1563
(
1998
).
80.
S. N.
Dorogovtsev
,
A. V.
Goltsev
, and
J. F. F.
Mendes
,
Rev. Mod. Phys.
80
,
1275
(
2008
).
81.
A.
Arenas
,
A.
Díaz-Guilera
,
J.
Kurths
,
Y.
Moreno
, and
C.
Zhou
,
Phys. Rep.
469
,
93
(
2008
).
82.
F.
Dörfler
and
F.
Bullo
,
Automatica
50
,
1539
(
2014
).
83.
S.
Ameli
,
F.
Shahbazi
,
M.
Karimian
, and
T.
Malakoutikhah
, e-print arXiv:1705.07875 (2017).
84.
S. H.
Strogatz
,
Physica D
143
,
1
(
2000
).
85.
D. A.
Wiley
,
S. H.
Strogatz
, and
M.
Girvan
,
Chaos
16
,
015103
(
2006
).
86.
D. M.
Abrams
and
S. H.
Strogatz
,
Phys. Rev. Lett.
93
,
174102
(
2004
).
87.
T.
Kotwal
,
X.
Jiang
, and
D. M.
Abrams
,
Phys. Rev. Lett.
119
,
264101
(
2017
).
88.
R.
Barber
,
M.
Horra
, and
J.
Crespo
,
IFAC Proc. Vol.
46
,
250
(
2013
).
89.
C. K.
Volos
,
I. M.
Kyprianidis
, and
I. N.
Stouboulos
,
Rob. Auton. Syst.
61
,
1314
(
2013
).
90.
F.
Candelas
,
G. J.
García
,
S.
Puente
,
J.
Pomares
,
C.
Jara
,
J.
Pérez
,
D.
Mira
, and
F.
Torres
,
IFAC-PapersOnLine
48
,
105
(
2015
).
91.
L.
de Castro
,
B.
Lago
, and
F.
Mondaini
,
J. Appl. Math. Phys.
3
,
631
(
2015
).
92.
A. D.
Pano-Azucena
,
J.
de Jesus Rangel-Magdaleno
,
E.
Tlelo-Cuautle
, and
A.
de Jesus Quintas-Valles
,
Nonlinear Dyn.
87
,
2203
(
2017
).
93.
A.
Kaboudia
,
E. M.
Cherry
, and
F. H.
Fenton
,
Sci. Adv.
5
,
aav6019
(
2019
).
94.
A.
Kaboudian
,
E. M.
Cherry
, and
F. H.
Fenton
,
Chaos Solitons Fractals
121
,
6
(
2019
).
95.
D. J.
Christini
,
M. L.
Riccio
,
C. A.
Culianu
,
J. J.
Fox
,
A.
Karma
, and
R. F.
Gilmour
, Jr.,
Phys. Rev. Lett.
96
,
104101
(
2006
).
96.
W.-J.
Rappel
,
F.
Fenton
, and
A.
Karma
,
Phys. Rev. Lett.
83
,
456
(
1999
).
97.
A.
Garzón
,
R. O.
Grigoriev
, and
F. H.
Fenton
,
Phys. Rev. E
80
,
021932
(
2009
).
98.
A.
Garzón
,
R. O.
Grigoriev
, and
F. H.
Fenton
,
Phys. Rev. E
84
,
041927
(
2011
).
99.
A.
Garzón
,
R. O.
Grigoriev
, and
F. H.
Fenton
,
Chaos
24
,
033124
(
2014
).
100.
E. M.
Cherry
,
Chaos
27
,
093902
(
2017
).
101.
See http://www.nimbios.org/workinggroups/WG_arrhythmias
for useful discussions on “Prediction and control of cardiac alternans
.”

Supplementary Material