Fireflies: New software for interactively exploring dynamical systems using GPU computing

In non-linear systems, where explicit analytic solutions usually can't be found, visualisation is a powerful approach which can give insights into the dynamical behaviour of models; it is also crucial for teaching this area of mathematics. In this paper we present new software, Fireflies, which exploits the power of graphical processing unit (GPU) computing to produce spectacular interactive visualisations of arbitrary systems of ordinary differential equations. In contrast to typical phase portraits, Fireflies draws the current position of trajectories (projected onto 2D or 3D space) as single points of light, which move as the system is simulated. Due to the massively parallel nature of GPU hardware, Fireflies is able to simulate millions of trajectories in parallel (even on standard desktop computer hardware), producing"swarms"of particles that move around the screen in real-time according to the equations of the system. Particles that move forwards in time reveal stable attractors (e.g. fixed points and limit cycles), while the option of integrating another group of trajectories backwards in time can reveal unstable objects (repellers). Fireflies allows the user to change the parameters of the system as it is running, in order to see the effect that they have on the dynamics and to observe bifurcations. We demonstrate the capabilities of the software with three examples: a two-dimensional"mean field"model of neuronal activity, the classical Lorenz system, and a 15-dimensional model of three interacting biologically realistic neurons.


Introduction
Many mathematical models are described by non-linear ordinary differential equations (ODEs).It is therefore important to be able to visualise solutions to these equations, and study how their dynamic behaviour changes when parameter values are modified.The theory of bifurcation analysis provides a rigorous mathematical framework for understanding how the stability of fixed points and limit cycles change as the parameters of the system are varied.Bifurcation theory can be combined with numerical continuation in order to trace the boundaries in parameter space that separate different dynamical regimes, and several software packages exist for this purpose, such as AUTO [1] and MATCONT [2].Before applying numerical continuation tools, in many cases it is useful to first get a rough understanding of how the phase space and system behaviour for particular parameter values, since the initial points for continuation (such as approximate fixed point and limit cycle positions) must be known.This is where more qualitative investigation (such as visualisation of the system) can be useful.
In this paper we introduce new software, "Fireflies" 1 for the dynamic visualisation of ODE solutions.Instead of the traditional method of showing complete trajectories in phase space, Fireflies presents the user with a twoor three-dimensional view of a cloud of moving particles.Each particle represents the position in state space of one trajectory at the current point in time, and as the simulation runs the particles move around the screen ac-cording to the equations of the system.Different coloured groups of particles can be given different ranges of initial conditions, and each group can be either integrated forwards or backwards in time.Forwards moving particles illuminate attractors (e.g.stable fixed points and limit cycles) of the system, while backwards moving ones illuminate repelling (unstable) objects.By watching the motion of the particles, it is also possible to see features such as saddle separatrices, and how the speed of movement varies around limit cycles.Finally, Fireflies can also generate "dynamic" bifurcation diagrams, in which one of the visualised dimensions corresponds to a chosen parameter and each particle is given a different value for this parameter.
When only the current position of each trajectory is shown, a very large number of particles must be used in order for the structure of the system to be visible -we find that several million are required for good results.Numerically integrating this many equations quickly enough to display an interactive animation is not typically possible using the limited parallelism of traditional central processing units (CPU), but is an ideal problem for so-called "massively parallel" graphical processing units (GPUs).Unlike CPUs, which typically contain a small number (up to 18, at time of writing) of largely independent processing cores, GPUs consist of hundreds or thousands of cores, all of which execute the same code at the same time.Nowadays GPUs, which were originally developed to accelerate rendering in 3D graphics applications, sit alongside the central processing unit (CPU) in almost all modern computers, including tablets and smartphones.Despite their origins in the video games industry, the use of GPUs is rapidly becoming an essential part of many scientific computing applications, due to their low cost, widespread availability, and potential for large efficiency gains in certain types of computations.Fireflies consists of a user friendly graphical interface which can be used to produce simulations of N-dimensional systems of ODEs.Thanks to the power of GPU computing, these simulations can contain millions of particles while still running quickly enough to be interactive, allowing the user to change parameter values and immediately observe the effect of this on the particles' motion.
In this paper we will demonstrate the capabilities of Fireflies by showing visualisations of three example systems of ODEs.We will show how using Fireflies to study these systems can give insights into their behaviour and how their dynamics depends on their parameters.Note that the figures show static screenshots from moving simulations; in order to properly appreciate the capabilities of the software we recommend viewing the movies included in the Supplementary Material.

Example Systems 2.1 A Two-Dimensional Model of the Basal Ganglia
In this section we show an example of a two dimensional system of nonlinear ODEs, corresponding to a simple model of neuronal activity.We show how Fireflies makes the dynamics of this system clear, and how the bifurcations of the system can be observed in an exciting new way.Although the system is two-dimensional, we also show how Fireflies can be used to create an interactive three-dimensional bifurcation diagram, by extending the visualisation into parameter space.
Based on earlier modelling work [3,4] we developed a mathematical model of activity in two connected regions in the brain: the subthalamic nucleus (STN) and external globus pallidus (GPe) [5].These regions are both part of the "basal ganglia", a group of nuclei that are known to be involved in movement control and that are severely affected by Parkinson's Disease.A current subject of much interest is the fact that the Parkinsonian basal ganglia show a much larger degree of rhythmically modulated (oscillatory) activity [6].We therefore use our model to study the conditions under which oscillations can appear in the STN-GPe network.
The equations of our model are: In these equations, which are based on the Wilson-Cowan formulation [7], the variables x and y correspond to the average level of neuronal activity in the STN and GPe populations, respectively.The populations are coupled to each other and themselves through chemical connections (synapses), and the strength of these connections are given by the parameters w gg , w gs , w ss and w sg (so, for example, w sg is the strength of the connection from STN to GPe).The STN population is excitatory, meaning it acts to increase the activity of its synaptic targets, while the GPe is inhibitory, meaning it acts to decrease the activity of its synaptic targets.The functions Z s and Z g are monotonically increasing sigmoid curves that describe how each population responds to synaptic input, and the parameters τ s and τ g are time constants that determine how quickly the activity in each population changes.
Here we will show how Fireflies can be used to investigate how the behaviour of the system changes with one of its parameters.The parameter we will vary is w ss , which corresponds to the strength of self-excitation in the STN.[3] demonstrated that if the neurons within the STN are able to excite each other (i.e.w ss > 0) then the STN-GPe circuit is able to generate oscillations, although currently there is no known biological mechanism for such self-excitation.The simulations shown all contain two "particle groups", one of which is integrated forwards in time and other backwards.Green and pink particles are used to render the forwards and backwards particle groups, respectively.Currently Fireflies supports only one solver method, 4 th order Runge-Kutta, with a constant time step that can be adjusted by the user during the simulation.All particles are given random initial conditions, with values of the system variables drawn from a uniform random distribution on the interval (0, 1); the values 0 and 1 correspond to the minimum and maximum levels of population activity respectively.Any trajectories that leave this square region, or which have not been reset for more than time T max , are reset to a new random set of initial conditions.
Figure 1A shows a screenshot of the simulation when w ss = 0, corresponding to the situation with no STN self-excitation.Under these conditions there is a single stable fixed point, which all of the green (forwards time) particles spiral in towards.The pink (backwards time) particles are difficult to see in this image, as they move very quickly out of the bounds of the system (to infinity) and are constantly being reset to new random initial conditions.It is interesting to note from this image that the density of particles approaching the spiral is not uniform, and has several discontinuities where the density suddenly changes.We suspected that this phenomena is a result of the boundaries of the phase space, and confirmed this with a traditional phase portrait where all the trajectories began on the edges of phase space (Figure 1B).The borders where particle density changes quickly correspond to the trajectories that begin in the four corners of phase space (red lines in Figure 1B); these four trajectories divide the phase plane into regions that "funnel" trajectories into the spiral.
The slider controls (visible in the bottom left of Figure 1A) can be used to explore the effects of changing the parameters of the system.When the value of a parameter is changed using its slider, the dynamics of the system change immediately and the difference can be clearly seen in the movement of the particles.Figure 1C shows the STN/GPe system after the strength of STN self-excitation (w ss ) has been increased to 4.9, taking the system through a fold of limit cycles bifurcation.For this parameter value, a pair of limit cycles (stable and unstable) now encircle the original stable fixed point.This illustrates the purpose of the group of particles that are integrated backwards (pink): these particles are attracted to unstable objects, here they reveal the location of the unstable limit cycle.As w ss is increased further, the unstable limit cycle shrinks around the fixed point and eventually disappears in a subcritical Andronov-Hopf bifurcation, at which point the fixed point becomes unstable.After this bifurcation the limit cycle is globally stable, but as w ss is increased further the period of oscillation increases.This can be seen on the visualisation in the form of particles "bunching up" and moving very slowly around one part of the cycle (Figure 1C).Eventually, a saddlenode on invariant circle (SNIC) bifurcation occurs, after which all trajectories approach a new, globally stable, fixed point.
Fireflies can also be used to generate three-dimensional animated bifurcation diagrams of two-dimensional systems.To see this, we redefine the system so that the bifurcation parameter (w ss ) is a new state variable, subject to dw ss /dt = 0.The initial condition range for this new state variable is set to be the range of parameter values that are of interest.With this configuration, each particle has its own value for w ss , and the space consists of a "continuum" of phase portraits, beautifully illustrating the dynamics of the system (Figure 2).The user can further investigate the effects of the other parameters on the system's dynamics by varying them interactively and observing how this changes the 3D bifurcation diagram.

The Lorenz Equations
In this example we use Fireflies to visualise the dynamics of the threedimensional Lorenz system.We describe in detail the results of various visualisations obtained for different values of a parameter of the system, and then produce an animated bifurcation diagram that summarises these results.
The classical Lorenz system consists of three coupled non-linear differential equations: ẋ = σ(y − x) Where σ, r and β are parameters of the system (σ, r, β > 0).These equations were originally studied as a simplified model of convection in the atmosphere [8], and a straightforward description of the different dynamical regimes that they can produce, including chaos, can be found in any textbook on nonlinear dynamics (e.g.[9, pp.311-320]).In this section we will demonstrate the ability of Fireflies to visualise three dimensional systems by exploring the behaviour of the Lorenz equations in the case where σ = 10 and β = 8  3 , while the parameter r is gradually increased from zero.For values of r less than 1, the origin is the only stable fixed point.Figure 3A shows Fireflies when r = 0.5, shortly after the simulation has started and all the particles are quickly approaching the origin.As was also seen in the 2D system described above, Fireflies shows how the space occupied by the cube of initial conditions is deformed around the fixed point as the particles move.
Note that we do not show any backwards-moving particles in this section, as they are not generally useful in visualisations of the Lorenz equations.This is because any unstable fixed points or cycles are of saddle type, which means that neither forward nor backward moving particles tend to them as t → ∞.
At r = 1 a supercritical pitchfork bifurcation occurs, and two new stable fixed points emerge from the origin.Figure 3B shows a screenshot from Fireflies with r = 5.With the new parameter value the origin becomes a saddle, and particles begin to move away from it in one of two directions, spiralling in to one of the new stable fixed points that were created in the bifurcation.To clearly show how trajectories spiral in to the fixed point, this simulation has two groups of particles (green and blue) with initial conditions that are drawn from separate small volumes in state space, both very close to the saddle at the origin, but on opposite sides of its incoming separatrix.The paths that the particles in these two groups take away from the saddle are very close to the saddle's outgoing separatrices.Trajectories starting at other points in state space all approach the spiral that is on the same side of the saddle's incoming separatrix as their initial position.
At r ≈ 13.926 a homoclinic bifurcation occurs, at which point the saddle's outgoing separatrices join up with its incoming ones, resulting in the creation of a pair of unstable limit cycles.For values of r beyond the bifurcation the separatrices have "crossed over", and trajectories that begin near the saddle make one cycle around their nearest spiral before returning back towards the saddle, and then spiralling in to the stable fixed point on the opposite side of the incoming separatrix, as shown in Figure 3C.Trajectories that start a bit further away from the saddle, however, rotate around their closest spiral at a much greater distance.When this rotation brings them close to the saddle's incoming separatrix, some of them split off and begin rotating around the opposite spiral.This unpredictable swapping, which corresponds to transient chaos, can happen many times before a particle finally spirals fully in to one of the two stable fixed points.Figure 3D shows one group of particles which all start closer to the top spiral than the bottom one, but a little further away from the origin than those in Figure 3C.Although most spiral in to the top fixed point, with each rotation a mass of particles splits off and switches to rotate around the bottom spiral.
Finally, at r ≈ 24.06 a strange attractor appears.Now, even though the two spirals remain stable until r ≈ 24.74, almost all trajectories flip backwards and forwards between the two spirals infinitely and chaotically.In Figure 3E, the two groups of particles that start near to the origin (green and blue) initially behave as in 3C, performing one large loop around their closest spiral before approaching quite close to the opposite spiral.The spirals are only very weakly stable now, however, and the trajectories spiral slowly away from them, instead of into them.After some time (as is just beginning to happen in 3E), the particles come far enough away from the spiral that some of them switch sides, beginning an endless series of such seemly random side swappings on the strange attractor.The figure also contains a third set of particles (bronze coloured) which start from a much wider set of initial conditions; these show the general shape of the strange attractor's surface.
By applying the same technique as in the previous section, we can also use Fireflies to generate an animated bifurcation diagram for the Lorenz equations.To achieve this, we make parameter r a state variable with ṙ = 0 and set up a two dimensional projection with axes r, y.Note that a three dimensional projection could also be used, but due to the perspective transformation the results are not shown as clearly in this case.Figure 4 shows a screenshot from a simulation using this configuration, where each particle's initial value of r is chosen from the interval (0 • • • 110).The supercritical pitchfork bifurcation is clearly visible at r = 1, and the appearance of a strange attractor at r ≈ 13.9 is marked by particles beginning to form a cloud that resembles noise.This cloud clearly has some very detailed structure, however, as can be seen more and more clearly as r increases.Several small windows of parameter values where the dynamics become regular are also visible, for example at r = 92.Again, it should be kept in mind that Figure 4 is a static screenshot of a moving animation of 8 million particles.The other two parameters of the system can be changed interactively, and the effect of this variation on the bifurcation diagram is seen immediately.

Coupled Hodgkin-Huxley Neurons
In this section we demonstrate the use of Fireflies to visualise a system that is considerably more complex than those presented above, consisting of 15 coupled ODEs that describe the activity in three synaptically connected neurons.
In 1963 Alan Hodgkin and Andrew Huxley received a Nobel prize for their discovery of the mechanism that allows neurons to generate action potentials (the electrical impulses, or "spikes", that neurons use to communicate with each other).More than 50 years later the general structure of the equations that Hodgkin and Huxley used to describe this mechanism still forms the basis of an enormous number of biologically realistic computational models of neuronal activity.Here, however, we will use the original equations and parameters, which specifically relate to biophysical activity in the giant axon of the squid [10].In this model, the electrical potential across a part of the neuron's membrane evolves according to the flow of sodium (Na) and potassium (K) ions through the membrane.The permeability of the membrane to these ions is controlled by gates, which open and close according to the membrane potential.Equations 6-9 describe the "classical" Hodgkin-Huxley model of neuronal dynamics in a population of N neurons.
Here V i is the membrane potential of the i th neuron and h i , m i and n i represent the average state of the gates on its ion channels (0 ≤ h i , m i , n i ≤ 1).The parameters C, ḡlk , ḡna , ḡk , e lk , e na and e k , along with the functions α X (.) and β X (.) (X ∈ {h, m, n}), are described in [10] or any computational neuroscience textbook.All parameters mentioned so far take the values given in [10].Parameter I i controls how much external current is injected into the cell -increasing this parameter causes a transition from quiescence to single action potential firing to repetitive firing.Finally, I i syn represents the total synaptic current that arises as a result of inputs from the other neurons in the model.If I i = I i syn = 0 then the membrane potential of neuron i approaches a fixed point at the "resting" potential of 0 (mV).In this example we use a network where the neurons are arranged in a loop, with each receiving synaptic input from only the previous neuron: The parameter ḡsyn is the maximum conductance (mS) of a synaptic connection (its "strength"); each synaptic connection has the same strength in this model.e syn is the synaptic equilibrium potential: if e syn > 0(mV ) then synaptic input is excitatory and acts to raise the membrane potential of the post-synaptic neuron above the resting state; in this example we set e syn = 10.For each neuron we consider a variable s i , where 0 < s i < 1, which acts as an indicator for when the neuron is spiking.This variable increases to 1 very rapidly (with time constant τ r ) when the neuron is firing an action potential, and then decays slowly (time constant τ d ) after a spike, as described by equation 11.Parameters σ and θ are the slope and shift respectively of the sigmoid function which is used to increase s i .
We begin by briefly showing a visualisation of the four dimensional phase space of a single independent neuron under varying strengths of current injection.For low values of I 1 all trajectories approach the fixed point which corresponds to the resting state.As I 1 is increased past I 1 ≈ 6.25 the resting state undergoes an Andronov-Hopf bifurcation and a stable limit cycle appears.This limit cycle corresponds to periodic spiking with a frequency that increases with I 1 .Figure 5 shows a screenshot of a simulation with a single neuron and I 1 = 10.The three dimensions of the projection are the three gating variables: h 1 , m 1 and n 1 .All the particles in this simulation eventually approach the stable limit cycle.However, there is also a spiralshaped cloud of particles (marked with an asterisk), which is made up of trajectories that start close to the formerly stable fixed point.For parameter value I 1 = 10 this fixed point is only weakly unstable, so trajectories nearby spend a long time oscillating with a low (gradually increasing) amplitude before they approach the limit cycle.
Next, we consider a ring of N = 3 neurons (Figure 6A), where each neuron receives an identical current injection that is strong enough to cause repetitive firing (I j = 10, j = 1, 2, 3).Such networks of coupled regular spiking neurons are typically able to synchronise their firing in either in-phase or anti-phase, depending on the strength and nature of synaptic connections [11].We set e syn = 10(mV ) and ḡsyn = 0.5(mS) to simulate the case of weak excitatory coupling.From running a few simulations of the system from random initial conditions in XPPAUT [12] it was clear that with these parameter values there was a very stable synchronous state, where all three neurons fired regularly with the same frequency and phase shift.Figure 6B shows a screenshot of Fireflies after this system has been simulated for long enough that all particles have settled down onto stable attractors.Each particle's position projected onto the space (V 1 , V 2 , V 3 ).The colour of each particle varies with its position,using a projection from the three dimensional co-ordinates to an RGB colour code.The prominenet solid white diagonal line is the synchronous limit cycle; most particles are attracted to this and continually move along the diagonal as the three neurons spike in unison.Figure 6C(i) shows the synchronous state as produced by XPPAUT.
In addition to the synchronous state, however, the visualisation also reveals four other stable limit cycles, three of which are easily visible with 500,000 particles.The three limit cycles appear as six coloured prongs in Figure 6B.By interactively moving around 3D space we were able to see that the particles were organised into three closed orbits, each of which included two differently coloured prongs.Each of these closed orbits corresponds to a limit cycle where two of the neurons are spiking and the third is silent.Simulations in XPPAUT revealed that each of these limit cycles corresponds to the state where the pre-synaptic neuron fires first, followed by the post-synaptic neuron, followed by a pause before the cycle repeats; this is shown in Figure 6C(ii).This detail about the order in which the neurons fire is not immediately visible in the Fireflies visualisation.
The other stable limit cycle that we have found has a relatively small basin of attraction, meaning that not many of the 500,000 particles approach it.It is also quite difficult to see in still screenshots, although some particles on this cycle are marked with an asterisk in Figure 6B.By following the path of these particles around the visualisation space, we can see that their trajectories correspond approach a stable limit cycle corresponding to the state where all three neurons spike in sequence.Interestingly, this order of this sequence is opposite to that of the direction of synaptic coupling (i.e. the firing sequence is neuron 1, neuron 3, neuron 2, ...).In order to verify that this was indeed a real limit cycle and not a numerical error in Fireflies, we wrote a short script to repeatedly (and sequentially) run simulations from random initial conditions in XPPAUT.The script recorded any sets of initial conditions that lead to solutions where all the neurons fired repetitively but non-synchronously.After performing a very large number of simulations this script had found only two sets of initial conditions that lead to the limit cycle that we had identified in Fireflies, and the resulting plots of V (t) are shown in Figure 6C(iii).
The example in this section demonstrates a case where Fireflies could reveal several stable dynamical regimes that were not immediately obvious from running individual simulations of the system from random initial conditions.Due to the fact that the system under study is composed of three identical elements, it is not surprising to find that multiple regimes corresponding to different symmetries are possible.Furthermore, it is likely that other limit cycles exist -for example there may be one in which all three neurons fire in sequence in order of synaptic coupling, rather than in the opposite order as we observed here.However, because none of the hundreds of thousands of particles in our simulations approached such attractors, it is likely that they are either unstable or have extremely small basins of attraction.

Discussion
To facilitate the qualitative study of systems of ODEs, software such as XP-PAUT [12], often includes the ability to plot trajectories in phase space that start from multiple sets of initial conditions, chosen either from a regular grid in phase space or stochastically.However, in systems with high dimensional phase space (many ODEs) or with finely structured dynamical regimes, a very large number of initial conditions must be used in order to have a high level of confidence that all stable attractors have been found.This can make the qualitative investigation a very slow process, especially if one wishes to also examine how the picture varies with parameter values.It is also very difficult to visualise many different trajectories on the same phase portrait, as large numbers of them will completely fill the phase space.We believe that Fireflies offers a useful and exciting new way to investigate dynamical systems.This approach is more intuitive than traditional phase portraits, as it presents dynamical systems as they are: dynamic.
In this paper we have demonstrated, through several examples, the power of Fireflies for qualitatively exploring the behaviour of dynamical systems.The visualisation of systems in explorable 2D or 3D space provide a new perspective on such systems that can help with intuitively understanding them.In addition to this, since Fireflies can simulate millions of trajectories in parallel, stable attractors can be discovered that might not have been found with other methods of qualitative investigation -as in section 3 (although clearly as the dimensionality of systems increases, the number of particles required to uniformly sample phase space with a given resolution increases exponentially).
Scientific computing is increasingly taking advantage of the power and availability of GPUs, and a range of techniques have emerged for using GPUs to solve problems in domains such as computational neuroscience [13] and fluid dynamics [14].To demonstrate the necessity of "massively parallel" GPU computing for producing Fireflies' visualisations, we compared the average time taken to compute a single integration step using a traditional CPU (single core and estimated quad core) to the time taken when integra- tion is performed on the GPU.Table 1 shows the results of this comparison.For interactive graphical applications, it is normally considered necessary to render at least 30 frames per second in order to give the user a smooth experience, which means that each frame must be generated in 33ms or less.Since the times shown in table 1 only represent the time taken to advance the simulation and do not include additional time needed to render particles to the screen, it would seem to be very difficult or impossible to run visualisations shown in this paper at an acceptable frame rate without using the GPU.Additionally, running the simulation on the CPU would also incur the significant additional overhead of passing the particles' current positions to the computer's graphical hardware for rendering; by performing integration on the GPU Fireflies minimizes the number of CPU-GPU memory transfers.
We have found that our new visualisation technique can be particularly useful in a teaching context, as it very easily demonstrates the behaviour that a particular set of equations produces.For example, visualisation of simple systems, such as the two-dimensional model of neuronal activity illustrated in section 2.1, can be used to show different bifurcations in a very intuitive way.For example, the parameter w ss is increased, Fireflies shows how particles move more and more slowly around one part of the limit cycle, becoming increasingly bunched together until finally being attracted into a new stable fixed point that appears in a SNIC bifurcation.
The ability of Fireflies to quickly simulate millions of trajectories at the same time is due to the fact that almost all processing occurs exclusively within the GPU, but this approach carries a number of disadvantages.In particular, due to the time taken to transfer data between GPU and CPU it is not possible to permanently record trajectories (e.g. for further analysis), however in future we plan to add the ability to record either a subset of trajectories, or all trajectories but with a recording interval that is greater than the time step.It should be noted, however, that this is not the main intended purpose of Fireflies, and other libraries for exploiting GPUs to solve ODEs, such as the Odeint library for C++ [15], may be more appropriate.Another limitation that arises from the use of the GPU is that Fireflies does not allow the equations of the system to contain branching (e.g.changing variable values in response to threshold crossing).This is because GPUs are based on a "single instruction multiple data" (SIMD) architecture, whereby each processing core executes the same instruction at the same time; due to this architecture branching causes a significant reduction in computation speed.
Fireflies is written in Python and utilises a number of cross-platform and open source libraries, namely NumPy [16], PyOpenGL2 , PyOpenCL [17], Py-Side 3 and Mako 4 .The software has been tested in Linux and Windows, and should also work in Apple OSX.The full source code for Fireflies, along with instructions for its use, can be found at https://bitbucket.org/rmerrison/fireflies.more Render Techniques and Particle Groups: • Each of the N State Variables corresponds to one of the system's ODEs and comprises a mathematical expression describing the variable's time derivative (right-hand side), along with the bounds of the variable.If a trajectory passes outside of bounds of any state variable during simulation it is reset to new random initial conditions.
• For each Parameter the user specifies a default value and the range of allowable values for the parameter.
• A Render Technique species a method for transforming the current position of each trajectory in phase space into a position of a particle on the screen.At present, particles can either be drawn using a simple 2D projection or a perspective 3D projection, and when specifying a projection the user must select which state variable corresponds to each of the (two or three) projection dimensions.The colour of rendered particles is also specified as part of a render technique; this can either be manually specified or varied automatically based on the position of the particle in 2D or 3D space.
• Each Particle Group consists of P particles that share the same render technique and range of initial conditions.The user specifies: the number of particles; the render technique that should be used to draw them; whether the particles in the group should be simulated forwards or backwards in time; and, for each state variable, a range of valid initial conditions.These initial condition ranges define a cube in phase space from which a position for each particle in the group is chosen when the particle is first initialized (or reset as a result of going out of bounds).
These objects that describe a system of equations are specified using the software's graphical user interface (Figure S1).Once this has been done, Fireflies initializes the system before entering its main execution loop, which consists of repeatedly updating the particles' positions and rendering them to the screen.An overview of how the software is organised is shown in Figure S2.The following sections describe each stage of processing in more detail.

A.1 Initialization
When a system is first created, or subsequently modified, the GPU must be configured to simulate it.This involves populating the GPU's memory with the particles' initial positions and compiling and uploading programs to update and render the system onto the GPU.Setting up the particles' initial positions is straightforward: if the system is N -dimensional and consists of a total of P particles, then an array with space for N • P single-precision floating point values is created and filled with random initial conditions (subject to the initial condition ranges specified by each particle's Group) -this array is then transferred into the GPU's memory.Note that there are two possible ways of organising this memory: the first involves storing the N values associated with the first particle, followed by those for the second particle, and so on (we call this "row-major" order), while the second method involves storing the P values corresponding to the particles' positions in the first dimension, followed by the values for the second dimension, etc (we call this "column-major" order).At present Fireflies always uses row-major order, however it is likely that due to the way in which GPU memory access operates, column-major may gives significant performance benefits; this is something we plan to investigate in a future version of the software.
In order to simulate the system on the GPU, a program (or "kernel") must be compiled that updates the position of a single particle.The general form of this kernel is independent of the particular system of equations being integrated, and involves reading the particle's current position from memory, determining its new position using a 4th order Runge-Kutta routine, and then writing the new position back into memory.The only part of the kernel that changes for different systems of equations is that which calculates the time derivative for a given point in state space (the right hand side of each ODE).This function is generated automatically using the definitions for the state variables given by the user.In order to convert the specified system into kernel source code (in the "OpenCL C" language), Fireflies uses a template source code file which is processed by the Mako templating library5 to produce the final kernel source.This source code is then compiled and uploaded onto the GPU using OpenCL.
In addition to the kernel, three other programs (or "shaders") must be compiled onto the GPU which perform the transformations necessary to render the particles to the screen.These shaders form a pipeline that converts from particle positions in state space into a series of squares (pairs of triangles) in two-dimensional screen space that can drawn by the graphics hardware.The different stages of this pipeline are the vertex shader, the geometry shader, and the fragment shader.
• The vertex shader takes as input the two or three components (depending on whether a 2D or 3D projection is required) of the particle's current position that correspond to the X, Y and Z directions of the scene to be rendered and outputs a two or three dimensional vector containing the position of the point transformed to be relative to an imaginary camera.This transformation is a simple linear transformation obtained by multiplying the vertex position by the "Model-View" matrix.In the case of a two-dimensional view, the camera always points directly at the plane containing the particles but it can be translated around and zoomed in and out.With a three-dimensional view, the user can interactively adjust the camera's position and orientation in space in order to "fly around" the system.As the user moves around the system, the Model-View matrix is updated accordingly.
• The geometry shader takes as input the camera-relative particle position from the vertex shader and first transforms this into two-dimensional position, using the "View-Projection" matrix.If the particle is to be rendered as a 2D projection then this transformation is straightforward, while for a 3D projection a standard perspective transformation matrix is used (for details of this see any 3D graphics programming resource).After the particle position has been projected into 2D, the geometry shader outputs four 2D vertices that correspond to a square that is centred on the particle position and facing the camera.
• The fragment shader is the final stage of the shader pipeline.When the graphics hardware renders the squares generated by the geometry shader on to the screen, it determines the colour of each pixel within the square by calling the fragment shader.The fragment shader generates an intensity value based on the distance from the centre of the square, producing a smooth circular particle shape.The actual colour output by the fragment shader is either specified explcitly by the particle's render technique, or varies linearly as a function of the particle's position in state space.The colour of the pixel on the screen is set to its current colour plus the colour contributed by the particle being drawn 6 .

A.2 Running the Simulation
Once the system has been specified and the GPU initialized with the kernel, shaders and initial particle positions, the simulation can be started by making the step size non-zero (both positive and negative step sizes are possible, to simulate the system either forwards or backwards in time).When the simulation is running, the main simulation loop repeatedly instructs the GPU to alternately execute the kernel and then the shader pipeline, causing the particles to be updated and drawn to the screen.The user can see, in real time, how the particles are moving around the state space and over time the system's attractors become clear.The user can explore the state space by changing the position of the camera in 2D or 3D space (effectively changing the Model-View matrix in GPU memory) using the mouse and keyboard.The system's parameters can also be changed during simulation by dragging sliders on the user interface, which updates the location in GPU memory where the corresponding parameter value is stored.
Some operations, such as checking for and resetting particles that have moved outside the bounds of the system to new random positions are not well-suited for GPU processing; Fireflies performs these operations on the CPU.Since copying data between the CPU and GPU is relatively slow, it is not practical to copy the positions of all of the particles back and forth at every time step.Instead, only a fraction of the particle position data is copied and processed by the CPU each frame.The CPU processes this batch of particles by checking if any of them have left the bounds of the system and chooses new initial conditions where necessary, before copying the affected particles' positions back to the GPU.This method ensures that simulations run at a smooth frame rate, although it means that particles are typically reset several time steps after their reset conditions are first satisfied.The user can adjust the size of CPU processing batches while the simulation is running, in order to find the setting that maximizes performance.initial period of unsynchronized firing, the particles all eventually approach of several limit cycles corresponding to different spiking patterns (see main text).Figure 2: 3D bifurcation diagram of the STN/GPe system under variation of the parameter w ss .The simulation contains 8 million particles, divided into a forwards in time group (green) and a backwards in time group (pink).The three visible bifurcations are marked with asterisks, from left to right: fold of limit cycles, subcritical Andronov-Hopf, saddle node on invariant circle.Note how the particles on the stable limit cycle become increasingly "bunched" at the top as the cycle approaches the SNIC bifurcation.The asterisks and direction arrows were added to the image manually and were not generated by Fireflies.Here, a state variable has been selected from the list of objects on the left hand side, allowing the properties of that state variable to be set using the panel on the right hand side.

Figure 1 :
Figure 1: Exploring the two dimensional STN/GPe model in Fireflies under variation of w ss .(A)-(D)show screenshots from Fireflies for different values of w ss , with 700,000 particles, half of which are integrated forwards in time (green) and half of which are integrated backwards (pink).(A) Full screenshot of the Fireflies window with w ss = 0.All forward-moving particles approach the globally stable spiral (in this still picture it is very difficult to see the backwards particles, which diverge quickly to infinity).(B) w ss = 4.9: A stable spiral, unstable limit cycle, and stable limit cycle.(C) w ss = 7.8:An unstable spiral and a stable limit cycle.(D) w ss = 11: An unstable node, saddle (not visible), and stable node.(E) Phase portrait (not generated by Fireflies) of the STN/GPe system with w ss = 0.All trajectories start from the borders of the phase plane and spiral in to the fixed point.The trajectories that start in the four corners of the phase plane are shown in pink.

Figure 3 :
Figure 3: Exploring the Lorenz equations by varying r.All particles are integrated forwards in time, different colours use different initial condition ranges.Asterisks approximately indicate the origin.(A) r = 0.5: All particles approach the stable fixed point at the origin.The complex shape is due to the boundaries of the initial conditions (x ∈ (−10, 10), y ∈ (−30, 30), z ∈ (0, 50)) and the shape of the fixed point's manifolds.(B) r = 5: The origin is now a saddle and two new stable fixed points have emerged.Both groups of particles have initial conditions near to the origin, but on opposite sides of its incoming separatrix (green: x, y ∈ (0, 0.01), z ∈ (0, 0.01), blue: x, y ∈ (−0.01, 0), z ∈ (0, 0.01)).Trajectories approach one of the spirals, based on which side of the separatrix they start on.(C-D) r = 15: Trajectories now switch between the two spirals for some time, before settling down and approaching one of them (transient chaos).(C) has initial conditions near the saddle, as in (B); these trajectories loop once around one spiral before approaching the other.(D) has initial conditions further away from the saddle; these trajectories can switch sides repeatedly.(E) r = 25: A strange attractor.Trajectories starting near the saddle (green and blue) begin by approaching one of the spirals, before slowly spiralling away and swapping sides chaotically.Another group of particles (bronze) with initial conditions spanning a large area show the shape of the attractor.27

Figure 4 :
Figure 4: A "live" bifurcation diagram for the Lorenz equations.The diagram was created by setting up a simulation with 8 million particles, each with a random fixed value for r, projected onto the (r, y) plane.The supercritical pitchfork at r = 1 is clearly visible, as is the appearance of a strange attractor at r ≈ 13.9 and various transient windows of non-chaotic behaviour such as at r = 92.Axes added manually.

Figure 5 :
Figure 5: A screenshot of Fireflies running the single Hodgkin-Huxley neuron model with 500,000 particles.The position of each particle is projected into the 3D space (h 1 , m 1 , n 1 ).Particles are coloured according to their position on each axis.Asterisk indicates the approximate position of the (weakly) unstable fixed point, with a cloud of particles moving slowly away from it.Axes added manually.

Figure 6 :
Figure 6: A Hodgkin-Huxley model with three neurons spiking repetitively in response to current injection, connected by weak excitatory synapses.(A) An overview of the model's structure.(B) A screenshot of Fireflies running the model with 500,000 particles.The position of each particle is projected into the 3D space (V 1 , V 2 , V 3 ).Particles are coloured according to their position on each axis.Asterisk indicates a small group of particles that are close to the limit cycle where each neuron spikes in turn.Axes added manually.(C) Individual plots (from XPPAUT) of trajectories starting from three different sets of initial conditions, leading to different limit cycles.i: Synchronous spiking; ii: Two neurons spiking, one suppressed; iii: Neurons spiking in turn.

Figure
Figure S1: A screenshot of the full Fireflies window showing the panel where the system to be simulated (the Lorenz equations in this case) is defined.Here, a state variable has been selected from the list of objects on the left hand side, allowing the properties of that state variable to be set using the panel on the right hand side.

Figure
Figure S2: A schematic overview of Fireflies.The processes at the top of the figure occur during intialization, and involve compiling an OpenCL kernel to advance the simulation, and generating random initial conditions.The main (lower) part of the figure shows the program's main loop, during which particles are repeatedly moved and then rendered to the screen.Note that almost all of the processing in the main loop occurs entirely within the GPU, avoiding the overhead of transferring data between the GPU and CPU.The exception to this is the resetting procedure, which only processes a subset of particles each frame.

Table 1 :
The time taken to advance each of the systems shown in this paper by a single step (averaged over 1000 steps), when executed on a single core CPU and a GPU.The estimated best-case performance for a quad-core CPU was estimated by dividing the single core figure by four.The CPU used was an Intel Core i5-2500k (3.3Ghz) and the GPU was an Nvidia Geforce GTX 460.The C code for the CPU benchmark (available at www.bitbucket.org/rmerrison/fireflies) was adapted from the corresponding OpenCL kernel for each system, and was compiled using GCC 4.8.2 on Ubuntu 14.04 using the -Ofast option.