Modeling Fireflies in Sync

Background

If you are in a rural part of Tennessee or south-east Asia on a warm summer night, you might be treated to one of mother natures original light shows. Staring into the dark of night, you will see a flash of light from a firefly. It will soon be followed by a second, and then a third, and so on until the whole landscape is covered in these small specs of light. If you keep watching them you’ll start to see them go from random flashes to large patters of flashing. Ultimately they will all be flashing in synchronous. This is not a trick of the mind finding patters from nothing, but a mating habit of the fireflies.

A nice video showing this effect is here.

This is something that requires hundreds of different fireflies to work together without a defined leader. It turns out that we can model the entire phenomenon as a series of coupled equations using something called the Kuramoto Model.

Kuramoto Model

The Kuramoto model is a way to describe large numbers of coupled oscillators. You can read more about it here, but the core of the model is the following governing equation

$\frac{d\theta_i}{dt} = \omega_i+\frac{K}{N}\sum_{j=1}^N\sin(\theta_j - \theta_i)$

Where ω is the natural frequency of the oscillator, and K is the coupling coefficient.

With this one governing equation, we can simulate N different oscillators. Lets now see how the assumptions made in the Kuramoto oscillator match up with the system we want to simulate. The fireflies turn “on” and “off” at a regular cycle, oscillating between the two states. They are coupled through seeing each other light up and turn off. This only occurs in certain species of fireflies so they are biologically similar. The fact that there is no leader firefly matches up nicely with the decentralized nature of the model.

Matlab Implementation

I am going to go through two Matlab implementations. The first is to simulate the fireflies using the powerful tools available in Matlab, namely ODE45. The second implementation will attempt to simulate the fireflies using only basic tooling that would be available in C++ using only the cmath (math.h) library. This second simulation will be used to determine the system requirements we will need to implement the model on an Arduino and will focus on a forward Euler method.

ODE implementation

Implementing this model in Matlab is quite straightforward. Each firefly’s phase offset is a state in the model, and you just need to loop through all the different combinations of fireflies to get the updates for the states. This can be done with two nested for loops.

In the following graph we can see 10 “Fireflies” synchronize up with each other after each starting with random initial phase offset. While they do coalesce to two different equilibrium, both are exactly in phase with each other due to the periodic nature of the sine function.

We can now plot the signal from the fireflies at specific points below using the following equation

$y(t) = \sin(t+ \theta(t))$

Below we can see that the “Fireflies” all start randomly and come into sync. If we were to increase the coupling coefficient it would happen at a much quicker rate.

Forward Euler Implementation

Now that we have simulated the fireflies in Matlab, can we transition it over to a real world demonstration piece? I don’t think this process needs all the computing power of a system running Matlab so lets see if I can shrink it down to fit on a micro controller. The man impediment to this is replacing the use of ODE45 in solving the ordinary differential governing equation. Lets start by looking at a generic Taylor series

$f(a+\Delta t)=f(a)+\frac{f'(a)}{1!}\Delta t +\frac{f''(a)}{2!}\Delta t^2 +... +\frac{f^{n}(a)}{n!}\Delta t^n = \sum_ {i=1}^n \frac{f^i(a)}{i!}\Delta t^i$

This gives us a general n’th order way of calculating the update to the function. If we took the limit as n goes to infinity we would match the function exactly (at least for a continuous function). We have the exact first derivative of θ and could get the following derivatives using numerical methods but these derivatives would contain lots of high frequency noise. In order to avoid all that noise lets truncate the Taylor series after the first derivative.

$f(a+ \Delta t) = f(a) + f'(a)\Delta t$

Estimating Error

Lets define the error that is caused by this truncation of the Taylor series as E, which is the difference between our truncated version to the infinitely long version.

$E(a +\Delta t) = f(a+\Delta t)_{truncated} - f(a+\Delta t)_{truth}= f(a) + f'(a)\Delta t - \sum_{i=1}^\infty \frac{f^i(a)}{i!}\Delta t^i$

Because, by definition, the first two terms of the true value match the truncated value. These cancel out and we can write the error as

$E(a+\Delta t) = \sum_{i=3}^\infty \frac{f^i(a)}{i!}\Delta t^i$

If we ensure that our Δt is smaller than 1 we can see how the error is dominated by the initial terms in our error equation. This means that as we make the time step smaller, the error introduced by our truncation will asymptotically go to 0. This isn’t completely true as past a certain point the error begins increasing due to errors introduced by the finite precision used by the computer. In reality there is an “error well” where an “optimal-ish” step size will provide the smallest error but this is problem specific and often hard to determine in non trivial problems.

Forward Euler Matlab

This method of stepping forward a differential equation using the first derivative of it is called the first order forward Euler method. If we were to expand this out to more terms it would be an n’th order Euler method. Then if we break down our Δt into subsections and evaluate the derivative of the function at several of these different subsections we would arrive at the Runge-Kutta methods (RK). If we, through a clever use of our error, allow our algorithm to change the time step dynamically we get to adaptive time step RK methods, which is exactly what Matlab uses in its non-fixed time step ODE suite for functions like ODE45. Matlab’s algorithms are more advanced than just standard adaptive step RK methods but I wanted to show how the forward Euler equation is just a simpler version of ODE45 which we used above.

Numerical integration techniques are still an active research area. I wanted only to give an extremely shallow introduction to how they work. Note: you should never trust numerical integration blindly. If you want an example of how they can fail you please see last weeks post where I tried to use ODE45 on a nonlinear system.

Well that was a fun tangent and we now know that we can replace ODE45 with the following equation for suitably small Δt to step our theta forward in time.

$\theta(k+1)_i=\theta(k)_i + \dot{\theta}(k)_i\Delta t = \theta(k)_i + (\omega_i+\frac{K}{N}\sum_{j=1}^{N}\sin(\theta_j - \theta_i))\Delta t$

We can then vary Δt, the period, in the simulation and see how small it needs to be to converge. If we take the inverse of the Δt we get the frequency that the micro controller will need to be running at. Note: I assumed that the calculations take a small fraction of the period of the controller. If the math took a significant percentage of the period we would need to include it.

Physical Implementation

While simulations on Matlab are good to play around with different parameters; graphs can only be so engaging. In order to pique interest and engage people with the topic, I tried to turn this simulation into a demonstration piece. I implemented the forward Euler method written about above on an Arduino. The only extra step added is to light up an LED every time the “Fireflies” signal is positive. Below is a video of the demonstration piece.

While not as impressive of a demonstration as hoped, it was still an interesting project to whip together in an afternoon. All code used in this post can be found here.

Want more Gereshes?

If you want to receive the weekly Gereshes blog post directly to your email every Monday morning you can sign up for the newsletter here!

If you can’t wait for next weeks post and want some more Gereshes I suggest

Rollout of a rocket motor test stand

How to pump a swing using math

Guerrilla astronomy