Fireflies: New software for interactively exploring dynamical systems using GPU computing
FFireflies: New software for interactivelyexploring dynamical systems using GPUcomputing
Robert Merrison-Hort ∗ School of Computing, Electronics & Mathematics, PlymouthUniversity
Abstract
In non-linear systems, where explicit analytic solutions usuallycan’t be found, visualisation is a powerful approach which can giveinsights into the dynamical behaviour of models; it is also crucial forteaching this area of mathematics. In this paper we present new soft-ware, Fireflies, which exploits the power of graphical processing unit(GPU) computing to produce spectacular interactive visualisations ofarbitrary systems of ordinary differential equations. In contrast totypical phase portraits, Fireflies draws the current position of trajec-tories (projected onto 2D or 3D space) as single points of light, whichmove as the system is simulated. Due to the massively parallel natureof GPU hardware, Fireflies is able to simulate millions of trajectoriesin parallel (even on standard desktop computer hardware), producing“swarms” of particles that move around the screen in real-time accord-ing to the equations of the system. Particles that move forwards intime reveal stable attractors (e.g. fixed points and limit cycles), whilethe option of integrating another group of trajectories backwards in ∗ Email: [email protected] a r X i v : . [ c s . M S ] M a y ime can reveal unstable objects (repellers). Fireflies allows the userto change the parameters of the system as it is running, in order to seethe 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 clas-sical Lorenz system, and a 15-dimensional model of three interactingbiologically realistic neurons. Many mathematical models are described by non-linear ordinary differentialequations (ODEs). It is therefore important to be able to visualise solutionsto these equations, and study how their dynamic behaviour changes whenparameter values are modified. The theory of bifurcation analysis providesa rigorous mathematical framework for understanding how the stability offixed points and limit cycles change as the parameters of the system arevaried. Bifurcation theory can be combined with numerical continuationin order to trace the boundaries in parameter space that separate differentdynamical regimes, and several software packages exist for this purpose, suchas AUTO [1] and MATCONT [2]. Before applying numerical continuationtools, in many cases it is useful to first get a rough understanding of howthe phase space and system behaviour for particular parameter values, sincethe initial points for continuation (such as approximate fixed point and limitcycle 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” for the dynamic vi-sualisation of ODE solutions. Instead of the traditional method of showingcomplete trajectories in phase space, Fireflies presents the user with a two-or three-dimensional view of a cloud of moving particles. Each particle rep-resents the position in state space of one trajectory at the current point intime, and as the simulation runs the particles move around the screen ac- Named after the classic dynamical systems example of firefly phase synchronisation,and also because the visualisations produced (slightly) resemble swarms of fireflies.
In this section we show an example of a two dimensional system of non-linear ODEs, corresponding to a simple model of neuronal activity. Weshow how Fireflies makes the dynamics of this system clear, and how thebifurcations of the system can be observed in an exciting new way. Althoughthe system is two-dimensional, we also show how Fireflies can be used tocreate an interactive three-dimensional bifurcation diagram, by extendingthe visualisation into parameter space.Based on earlier modelling work [3, 4] we developed a mathematical modelof activity in two connected regions in the brain: the subthalamic nucleus(STN) and external globus pallidus (GPe) [5]. These regions are both partof the “basal ganglia”, a group of nuclei that are known to be involvedin movement control and that are severely affected by Parkinson’s Disease.A current subject of much interest is the fact that the Parkinsonian basalganglia show a much larger degree of rhythmically modulated (oscillatory)activity [6]. We therefore use our model to study the conditions under whichoscillations can appear in the STN-GPe network.The equations of our model are: τ s ˙ x = − x + Z s ( w ss x − w gs y + I ) (1) τ g ˙ y = − y + Z g ( − w gg y + w sg x ) (2)In these equations, which are based on the Wilson-Cowan formulation47], the variables x and y correspond to the average level of neuronal activityin the STN and GPe populations, respectively. The populations are coupledto each other and themselves through chemical connections (synapses), andthe 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 STNto GPe). The STN population is excitatory, meaning it acts to increase theactivity of its synaptic targets, while the GPe is inhibitory, meaning it actsto decrease the activity of its synaptic targets. The functions Z s and Z g aremonotonically increasing sigmoid curves that describe how each populationresponds to synaptic input, and the parameters τ s and τ g are time constantsthat determine how quickly the activity in each population changes.Here we will show how Fireflies can be used to investigate how the be-haviour of the system changes with one of its parameters. The parameterwe will vary is w ss , which corresponds to the strength of self-excitation inthe STN. [3] demonstrated that if the neurons within the STN are able toexcite each other (i.e. w ss >
0) then the STN-GPe circuit is able to generateoscillations, although currently there is no known biological mechanism forsuch self-excitation. The simulations shown all contain two “particle groups”,one of which is integrated forwards in time and other backwards. Green andpink particles are used to render the forwards and backwards particle groups,respectively. Currently Fireflies supports only one solver method, 4 th orderRunge-Kutta, with a constant time step that can be adjusted by the userduring the simulation. All particles are given random initial conditions, withvalues of the system variables drawn from a uniform random distribution onthe interval (0 , T max ,are reset to a new random set of initial conditions.Figure 1A shows a screenshot of the simulation when w ss = 0, correspond-ing to the situation with no STN self-excitation. Under these conditions thereis a single stable fixed point, which all of the green (forwards time) particlesspiral in towards. The pink (backwards time) particles are difficult to seein this image, as they move very quickly out of the bounds of the system5to infinity) and are constantly being reset to new random initial conditions.It is interesting to note from this image that the density of particles ap-proaching the spiral is not uniform, and has several discontinuities wherethe density suddenly changes. We suspected that this phenomena is a resultof the boundaries of the phase space, and confirmed this with a traditionalphase portrait where all the trajectories began on the edges of phase space(Figure 1B). The borders where particle density changes quickly correspondto the trajectories that begin in the four corners of phase space (red lines inFigure 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 usedto explore the effects of changing the parameters of the system. When thevalue of a parameter is changed using its slider, the dynamics of the systemchange immediately and the difference can be clearly seen in the movementof the particles. Figure 1C shows the STN/GPe system after the strengthof STN self-excitation ( w ss ) has been increased to 4.9, taking the systemthrough a fold of limit cycles bifurcation. For this parameter value, a pair oflimit cycles (stable and unstable) now encircle the original stable fixed point.This illustrates the purpose of the group of particles that are integratedbackwards (pink): these particles are attracted to unstable objects, here theyreveal the location of the unstable limit cycle. As w ss is increased further, theunstable limit cycle shrinks around the fixed point and eventually disappearsin a subcritical Andronov-Hopf bifurcation, at which point the fixed pointbecomes unstable. After this bifurcation the limit cycle is globally stable,but as w ss is increased further the period of oscillation increases. This can beseen on the visualisation in the form of particles “bunching up” and movingvery slowly around one part of the cycle (Figure 1C). Eventually, a saddle-node on invariant circle (SNIC) bifurcation occurs, after which all trajectoriesapproach a new, globally stable, fixed point.Fireflies can also be used to generate three-dimensional animated bifur-cation diagrams of two-dimensional systems. To see this, we redefine thesystem so that the bifurcation parameter ( w ss ) is a new state variable, sub-ject to dw ss /dt = 0. The initial condition range for this new state variable6s set to be the range of parameter values that are of interest. With thisconfiguration, each particle has its own value for w ss , and the space consistsof a “continuum” of phase portraits, beautifully illustrating the dynamics ofthe system (Figure 2). The user can further investigate the effects of theother parameters on the system’s dynamics by varying them interactivelyand observing how this changes the 3D bifurcation diagram. In this example we use Fireflies to visualise the dynamics of the three-dimensional Lorenz system. We describe in detail the results of variousvisualisations obtained for different values of a parameter of the system,and then produce an animated bifurcation diagram that summarises theseresults.The classical Lorenz system consists of three coupled non-linear differen-tial equations: ˙ x = σ ( y − x ) (3)˙ y = x ( r − z ) − y (4)˙ z = xy − βz (5)Where σ , r and β are parameters of the system ( σ, r, β > σ = 10and β = , while the parameter r is gradually increased from zero.For values of r less than 1, the origin is the only stable fixed point. Figure3A shows Fireflies when r = 0 .
5, shortly after the simulation has started andall the particles are quickly approaching the origin. As was also seen in the 2Dsystem described above, Fireflies shows how the space occupied by the cubeof initial conditions is deformed around the fixed point as the particles move.7ote that we do not show any backwards-moving particles in this section, asthey are not generally useful in visualisations of the Lorenz equations. Thisis because any unstable fixed points or cycles are of saddle type, which meansthat neither forward nor backward moving particles tend to them as t → ∞ .At r = 1 a supercritical pitchfork bifurcation occurs, and two new stablefixed points emerge from the origin. Figure 3B shows a screenshot fromFireflies with r = 5. With the new parameter value the origin becomes asaddle, 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 thebifurcation. To clearly show how trajectories spiral in to the fixed point, thissimulation has two groups of particles (green and blue) with initial conditionsthat are drawn from separate small volumes in state space, both very closeto 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 saddleare very close to the saddle’s outgoing separatrices. Trajectories starting atother points in state space all approach the spiral that is on the same side ofthe saddle’s incoming separatrix as their initial position.At r ≈ .
926 a homoclinic bifurcation occurs, at which point the saddle’soutgoing separatrices join up with its incoming ones, resulting in the creationof a pair of unstable limit cycles. For values of r beyond the bifurcation theseparatrices have “crossed over”, and trajectories that begin near the saddlemake one cycle around their nearest spiral before returning back towards thesaddle, and then spiralling in to the stable fixed point on the opposite side ofthe incoming separatrix, as shown in Figure 3C. Trajectories that start a bitfurther away from the saddle, however, rotate around their closest spiral at amuch greater distance. When this rotation brings them close to the saddle’sincoming separatrix, some of them split off and begin rotating around theopposite spiral. This unpredictable swapping, which corresponds to transientchaos, can happen many times before a particle finally spirals fully in to oneof the two stable fixed points. Figure 3D shows one group of particles whichall start closer to the top spiral than the bottom one, but a little furtheraway from the origin than those in Figure 3C. Although most spiral in to thetop fixed point, with each rotation a mass of particles splits off and switches8o rotate around the bottom spiral.Finally, at r ≈ .
06 a strange attractor appears. Now, even thoughthe two spirals remain stable until r ≈ .
74, almost all trajectories flipbackwards and forwards between the two spirals infinitely and chaotically.In Figure 3E, the two groups of particles that start near to the origin (greenand blue) initially behave as in 3C, performing one large loop around theirclosest spiral before approaching quite close to the opposite spiral. Thespirals are only very weakly stable now, however, and the trajectories spiralslowly away from them, instead of into them. After some time (as is justbeginning to happen in 3E), the particles come far enough away from thespiral that some of them switch sides, beginning an endless series of suchseemly random side swappings on the strange attractor. The figure alsocontains a third set of particles (bronze coloured) which start from a muchwider set of initial conditions; these show the general shape of the strangeattractor’s surface.By applying the same technique as in the previous section, we can also useFireflies to generate an animated bifurcation diagram for the Lorenz equa-tions. To achieve this, we make parameter r a state variable with ˙ r = 0and set up a two dimensional projection with axes r, y . Note that a threedimensional projection could also be used, but due to the perspective trans-formation the results are not shown as clearly in this case. Figure 4 shows ascreenshot from a simulation using this configuration, where each particle’sinitial value of r is chosen from the interval (0 · · · r = 1, and the appearance of astrange attractor at r ≈ . r increases. Severalsmall windows of parameter values where the dynamics become regular arealso visible, for example at r = 92. Again, it should be kept in mind thatFigure 4 is a static screenshot of a moving animation of 8 million particles.The other two parameters of the system can be changed interactively, andthe effect of this variation on the bifurcation diagram is seen immediately.9 Coupled Hodgkin-Huxley Neurons
In this section we demonstrate the use of Fireflies to visualise a system thatis considerably more complex than those presented above, consisting of 15coupled ODEs that describe the activity in three synaptically connected neu-rons.In 1963 Alan Hodgkin and Andrew Huxley received a Nobel prize for theirdiscovery of the mechanism that allows neurons to generate action potentials(the electrical impulses, or “spikes”, that neurons use to communicate witheach other). More than 50 years later the general structure of the equationsthat Hodgkin and Huxley used to describe this mechanism still forms thebasis of an enormous number of biologically realistic computational modelsof neuronal activity. Here, however, we will use the original equations andparameters, which specifically relate to biophysical activity in the giant axonof the squid [10]. In this model, the electrical potential across a part of theneuron’s membrane evolves according to the flow of sodium (Na) and potas-sium (K) ions through the membrane. The permeability of the membraneto these ions is controlled by gates, which open and close according to themembrane potential. Equations 6–9 describe the “classical” Hodgkin-Huxleymodel of neuronal dynamics in a population of N neurons. C ˙ V i = ¯ g lk ( e lk − V i ) + h i m i ¯ g na ( e na − V i ) + n i ¯ g k ( e k − V i ) + I isyn ( t ) + I i (6)˙ h i = α h ( V i )(1 − h i ) − β h ( V i ) h i (7)˙ m i = α m ( V i )(1 − m i ) − β m ( V i ) m i (8)˙ n i = α n ( V i )(1 − n i ) − β n ( V i ) n i (9) i = 1 , , ..., N 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 ≤ C , ¯ g lk , ¯ g na , ¯ g 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 computational10euroscience textbook. All parameters mentioned so far take the values givenin [10]. Parameter I i controls how much external current is injected into thecell - increasing this parameter causes a transition from quiescence to singleaction potential firing to repetitive firing. Finally, I isyn represents the totalsynaptic current that arises as a result of inputs from the other neurons in themodel. If I i = I isyn = 0 then the membrane potential of neuron i approachesa fixed point at the “resting” potential of 0 (mV). In this example we usea network where the neurons are arranged in a loop, with each receivingsynaptic input from only the previous neuron: I isyn ( t ) = ¯ g syn ( e syn − V i ) s ( i − modN (10)˙ s i = τ − r (1 + exp ( − σ ( V i − θ ))) − (1 − s i ) − τ − d s i (11)The parameter ¯ g syn is the maximum conductance (mS) of a synapticconnection (its “strength”); each synaptic connection has the same strengthin this model. e syn is the synaptic equilibrium potential: if e syn > mV )then synaptic input is excitatory and acts to raise the membrane potentialof 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 < τ r ) when the neuron is firingan action potential, and then decays slowly (time constant τ d ) after a spike,as described by equation 11. Parameters σ and θ are the slope and shiftrespectively of the sigmoid function which is used to increase s i .We begin by briefly showing a visualisation of the four dimensional phasespace of a single independent neuron under varying strengths of current in-jection. For low values of I all trajectories approach the fixed point whichcorresponds to the resting state. As I is increased past I ≈ .
25 the rest-ing state undergoes an Andronov-Hopf bifurcation and a stable limit cycleappears. This limit cycle corresponds to periodic spiking with a frequencythat increases with I . Figure 5 shows a screenshot of a simulation with asingle neuron and I = 10. The three dimensions of the projection are thethree gating variables: h , m and n . All the particles in this simulation11ventually approach the stable limit cycle. However, there is also a spiral-shaped cloud of particles (marked with an asterisk), which is made up oftrajectories that start close to the formerly stable fixed point. For parametervalue I = 10 this fixed point is only weakly unstable, so trajectories nearbyspend a long time oscillating with a low (gradually increasing) amplitudebefore they approach the limit cycle.Next, we consider a ring of N = 3 neurons (Figure 6A), where eachneuron receives an identical current injection that is strong enough to causerepetitive firing ( I j = 10 , j = 1 , , e syn = 10( mV ) and ¯ g syn = 0 . mS ) to simulate the case ofweak excitatory coupling. From running a few simulations of the systemfrom random initial conditions in XPPAUT [12] it was clear that with theseparameter values there was a very stable synchronous state, where all threeneurons fired regularly with the same frequency and phase shift. Figure 6Bshows a screenshot of Fireflies after this system has been simulated for longenough that all particles have settled down onto stable attractors. Eachparticle’s position projected onto the space ( V , V , V ). The colour of eachparticle varies with its position,using a projection from the three dimensionalco-ordinates to an RGB colour code. The prominenet solid white diagonalline is the synchronous limit cycle; most particles are attracted to this andcontinually 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 re-veals four other stable limit cycles, three of which are easily visible with500,000 particles. The three limit cycles appear as six coloured prongs inFigure 6B. By interactively moving around 3D space we were able to see thatthe particles were organised into three closed orbits, each of which includedtwo differently coloured prongs. Each of these closed orbits corresponds toa limit cycle where two of the neurons are spiking and the third is silent.Simulations in XPPAUT revealed that each of these limit cycles correspondsto the state where the pre-synaptic neuron fires first, followed by the post-12ynaptic neuron, followed by a pause before the cycle repeats; this is shownin Figure 6C(ii). This detail about the order in which the neurons fire is notimmediately visible in the Fireflies visualisation.The other stable limit cycle that we have found has a relatively small basinof attraction, meaning that not many of the 500,000 particles approach it.It is also quite difficult to see in still screenshots, although some particleson this cycle are marked with an asterisk in Figure 6B. By following thepath of these particles around the visualisation space, we can see that theirtrajectories correspond approach a stable limit cycle corresponding to thestate where all three neurons spike in sequence. Interestingly, this order ofthis 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 verifythat 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 fromrandom initial conditions in XPPAUT. The script recorded any sets of initialconditions that lead to solutions where all the neurons fired repetitively butnon-synchronously. After performing a very large number of simulations thisscript had found only two sets of initial conditions that lead to the limit cyclethat we had identified in Fireflies, and the resulting plots of V ( t ) are shownin Figure 6C(iii).The example in this section demonstrates a case where Fireflies could re-veal several stable dynamical regimes that were not immediately obvious fromrunning individual simulations of the system from random initial conditions.Due to the fact that the system under study is composed of three identicalelements, it is not surprising to find that multiple regimes corresponding todifferent symmetries are possible. Furthermore, it is likely that other limitcycles exist - for example there may be one in which all three neurons fire insequence in order of synaptic coupling, rather than in the opposite order aswe observed here. However, because none of the hundreds of thousands ofparticles in our simulations approached such attractors, it is likely that theyare either unstable or have extremely small basins of attraction.13 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 thatstart from multiple sets of initial conditions, chosen either from a regular gridin phase space or stochastically. However, in systems with high dimensionalphase space (many ODEs) or with finely structured dynamical regimes, avery large number of initial conditions must be used in order to have a highlevel of confidence that all stable attractors have been found. This can makethe qualitative investigation a very slow process, especially if one wishes toalso examine how the picture varies with parameter values. It is also verydifficult to visualise many different trajectories on the same phase portrait,as large numbers of them will completely fill the phase space. We believethat Fireflies offers a useful and exciting new way to investigate dynamicalsystems. This approach is more intuitive than traditional phase portraits, asit presents dynamical systems as they are: dynamic -Ofast option.tion is performed on the GPU. Table 1 shows the results of this comparison.For interactive graphical applications, it is normally considered necessary torender at least 30 frames per second in order to give the user a smooth ex-perience, 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 advancethe simulation and do not include additional time needed to render particlesto the screen, it would seem to be very difficult or impossible to run visual-isations shown in this paper at an acceptable frame rate without using theGPU. Additionally, running the simulation on the CPU would also incur thesignificant additional overhead of passing the particles’ current positions tothe computer’s graphical hardware for rendering; by performing integrationon the GPU Fireflies minimizes the number of CPU-GPU memory transfers.We have found that our new visualisation technique can be particularlyuseful in a teaching context, as it very easily demonstrates the behaviour thata particular set of equations produces. For example, visualisation of simplesystems, such as the two-dimensional model of neuronal activity illustrated insection 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 particlesmove more and more slowly around one part of the limit cycle, becomingincreasingly bunched together until finally being attracted into a new stable15xed point that appears in a SNIC bifurcation.The ability of Fireflies to quickly simulate millions of trajectories at thesame time is due to the fact that almost all processing occurs exclusivelywithin the GPU, but this approach carries a number of disadvantages. Inparticular, due to the time taken to transfer data between GPU and CPU itis 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 oftrajectories, or all trajectories but with a recording interval that is greaterthan the time step. It should be noted, however, that this is not the mainintended purpose of Fireflies, and other libraries for exploiting GPUs to solveODEs, 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 doesnot allow the equations of the system to contain branching (e.g. changingvariable values in response to threshold crossing). This is because GPUs arebased on a “single instruction multiple data” (SIMD) architecture, wherebyeach processing core executes the same instruction at the same time; dueto this architecture branching causes a significant reduction in computationspeed.Fireflies is written in Python and utilises a number of cross-platform andopen source libraries, namely NumPy [16], PyOpenGL , PyOpenCL [17], Py-Side and Mako . The software has been tested in Linux and Windows, andshould also work in Apple OSX. The full source code for Fireflies, along withinstructions for its use, can be found at https://bitbucket.org/rmerrison/fireflies. Acknowledgement
I am extremely grateful to Roman Borisyuk for the valuable input thathe gave during the development of Fireflies and the preparation of thismanuscript. http://pyopengl.sourceforge.net/ eferences [1] E. J. Doedel and B. E. Oldeman, “Auto-07p: Continuation and bifurca-tion software for ordinary differential equations,” tech. rep., ConcordiaUniversity, Montreal, 2012.[2] A. Dhooge, W. Govaerts, and Y. A. Kuznetsov, “Matcont: A matlabpackage for numerical bifurcation analysis of odes,” ACM Transactionson Mathematical Software , vol. 29, no. 2, pp. 141–164, 2003.[3] A. Gillies, D. Willshaw, and Z. Li, “Subthalamic-pallidal interactionsare critical in determining normal and abnormal functioning of the basalganglia,”
Proceedings of the Royal Society of London. Series B: Biolog-ical Sciences , vol. 269, no. 1491, pp. 545–551, 2002.[4] A. Pavlides, S. John Hogan, and R. Bogacz, “Improved conditions forthe generation of beta oscillations in the subthalamic nucleus-globuspallidus network,”
European Journal of Neuroscience , vol. 36, no. 2,pp. 2229–2239, 2012.[5] R. J. Merrison-Hort, N. Yousif, F. Njap, U. G. Hofmann, O. Burylko,and R. Borisyuk, “An interactive channel model of the basal ganglia:Bifurcation analysis under healthy and parkinsonian conditions,”
TheJournal of Mathematical Neuroscience , vol. 3, no. 1, 2013.[6] T. Boraud, P. Brown, J. Goldberg, A. Graybiel, and P. Magill, “Oscil-lations in the Basal Ganglia: The good, the bad, and the unexpected,”in
The Basal Ganglia VIII (J. Bolam, C. Ingham, and P. Magill, eds.),ch. 1, pp. 1–24, Springer Science and Business Media, 2005.[7] H. Wilson and J. Cowan, “Exictatory and inhibitory interactions inlocalized populations of model neurons,”
Biophysical Journal , vol. 12,pp. 1–24, 1972.[8] E. N. Lorenz, “Deterministic nonperiodic flow,”
Journal of the Atmo-spheric Sciences , vol. 20, pp. 130–141, 1963.179] S. H. Strogatz,
Nonlinear Dynamics and Chaos . Westview Press, 1994.[10] A. L. Hodgkin and A. F. Huxley, “A quantitative description of mem-brane current and its application to conduction and excitation in nerve,”
The Journal of Physiology , vol. 117, no. 4, p. 500, 1952.[11] R. E. Mirollo and S. H. Strogatz, “Synchronization of pulse-coupledbiological oscillators,”
SIAM Journal on Applied Mathematics , vol. 50,no. 6, pp. 1645–1662, 1990.[12] B. Ermentrout,
Simulating, Analyzing, and Animating Dynamical Sys-tems: A Guide to XPPAUT for Researchers and Students . Philadelpha,PA, USA: SIAM, 2002.[13] R. Brette and D. F. M. Goodman, “Simulating spiking neural networkson gpu,”
Network: Computation in Neural Systems , vol. 23, no. 4,pp. 167–182, 2012.[14] M. J. Harris,
GPU Gems , ch. Fast Fluid Dynamics Simulation on theGPU, pp. 637–665. Addison-Wesley, 2004.[15] K. Ahnert and M. Mulansky, “Odeint - solving ordinary differentialequations in c++.” http://arxiv.org/abs/1110.3397v1, 2011.[16] S. van der Walt, S. Colbert, and G. Varoquaux, “The numpy array:A structure for efficient numerical computation,”
Computing in ScienceEngineering , vol. 13, pp. 22–30, March 2011.[17] A. Kl¨ockner, N. Pinto, Y. Lee, B. Catanzaro, P. Ivanov, and A. Fasih,“PyCUDA and PyOpenCL: A Scripting-Based Approach to GPU Run-Time Code Generation,”
Parallel Computing , vol. 38, no. 3, pp. 157–174,2012.
Appendix A Implementation Details
In order to simulate a system of ODEs in Fireflies, the user must specifythe N State Variables , along with the
Parameters of the system, and one or18ore
Render Techniques and
Particle Groups : • Each of the N State Variables corresponds to one of the system’s ODEsand comprises a mathematical expression describing the variable’s timederivative (right-hand side), along with the bounds of the variable.If a trajectory passes outside of bounds of any state variable duringsimulation it is reset to new random initial conditions. • For each
Parameter the user specifies a default value and the range ofallowable values for the parameter. • A Render Technique species a method for transforming the currentposition of each trajectory in phase space into a position of a particleon the screen. At present, particles can either be drawn using a simple2D projection or a perspective 3D projection, and when specifying aprojection the user must select which state variable corresponds to eachof the (two or three) projection dimensions. The colour of renderedparticles is also specified as part of a render technique; this can eitherbe manually specified or varied automatically based on the position ofthe particle in 2D or 3D space. • Each
Particle Group consists of P particles that share the same ren-der technique and range of initial conditions. The user specifies: thenumber of particles; the render technique that should be used to drawthem; whether the particles in the group should be simulated forwardsor backwards in time; and, for each state variable, a range of valid ini-tial conditions. These initial condition ranges define a cube in phasespace from which a position for each particle in the group is chosenwhen the particle is first initialized (or reset as a result of going out ofbounds).These objects that describe a system of equations are specified using thesoftware’s graphical user interface (Figure S1). Once this has been done,Fireflies initializes the system before entering its main execution loop, whichconsists of repeatedly updating the particles’ positions and rendering them to19he screen. An overview of how the software is organised is shown in FigureS2. 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 beconfigured to simulate it. This involves populating the GPU’s memory withthe particles’ initial positions and compiling and uploading programs to up-date and render the system onto the GPU. Setting up the particles’ initialpositions is straightforward: if the system is N -dimensional and consists of atotal of P particles, then an array with space for N · P single-precision float-ing point values is created and filled with random initial conditions (subjectto the initial condition ranges specified by each particle’s Group ) - this arrayis then transferred into the GPU’s memory. Note that there are two possi-ble ways of organising this memory: the first involves storing the N valuesassociated with the first particle, followed by those for the second particle,and so on (we call this “row-major” order), while the second method in-volves storing the P values corresponding to the particles’ positions in thefirst dimension, followed by the values for the second dimension, etc (we callthis “column-major” order). At present Fireflies always uses row-major or-der, however it is likely that due to the way in which GPU memory accessoperates, column-major may gives significant performance benefits; this issomething 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 generalform of this kernel is independent of the particular system of equations beingintegrated, and involves reading the particle’s current position from memory,determining its new position using a 4th order Runge-Kutta routine, andthen writing the new position back into memory. The only part of the kernelthat changes for different systems of equations is that which calculates thetime derivative for a given point in state space (the right hand side of eachODE). This function is generated automatically using the definitions for thestate variables given by the user. In order to convert the specified system20nto kernel source code (in the “OpenCL C” language), Fireflies uses a tem-plate source code file which is processed by the Mako templating library to produce the final kernel source. This source code is then compiled anduploaded onto the GPU using OpenCL.In addition to the kernel, three other programs (or “shaders”) must becompiled onto the GPU which perform the transformations necessary to ren-der the particles to the screen. These shaders form a pipeline that convertsfrom particle positions in state space into a series of squares (pairs of trian-gles) in two-dimensional screen space that can drawn by the graphics hard-ware. The different stages of this pipeline are the vertex shader, the geometryshader, and the fragment shader. • The vertex shader takes as input the two or three components (de-pending on whether a 2D or 3D projection is required) of the particle’scurrent position that correspond to the X, Y and Z directions of thescene to be rendered and outputs a two or three dimensional vectorcontaining the position of the point transformed to be relative to animaginary camera. This transformation is a simple linear transforma-tion obtained by multiplying the vertex position by the “Model-View”matrix. In the case of a two-dimensional view, the camera always pointsdirectly at the plane containing the particles but it can be translatedaround and zoomed in and out. With a three-dimensional view, theuser can interactively adjust the camera’s position and orientation inspace in order to “fly around” the system. As the user moves aroundthe system, the Model-View matrix is updated accordingly. • The geometry shader takes as input the camera-relative particle posi-tion from the vertex shader and first transforms this into two-dimensionalposition, using the “View-Projection” matrix. If the particle is to berendered as a 2D projection then this transformation is straightfor-ward, while for a 3D projection a standard perspective transformationmatrix is used (for details of this see any 3D graphics programmingresource). After the particle position has been projected into 2D, the • The fragment shader is the final stage of the shader pipeline. Whenthe graphics hardware renders the squares generated by the geometryshader on to the screen, it determines the colour of each pixel withinthe square by calling the fragment shader. The fragment shader gen-erates an intensity value based on the distance from the centre of thesquare, producing a smooth circular particle shape. The actual colouroutput by the fragment shader is either specified explcitly by the par-ticle’s render technique, or varies linearly as a function of the particle’sposition in state space. The colour of the pixel on the screen is setto its current colour plus the colour contributed by the particle beingdrawn . 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 makingthe step size non-zero (both positive and negative step sizes are possible, tosimulate the system either forwards or backwards in time).When the simulation is running, the main simulation loop repeatedlyinstructs the GPU to alternately execute the kernel and then the shaderpipeline, causing the particles to be updated and drawn to the screen. Theuser can see, in real time, how the particles are moving around the statespace and over time the system’s attractors become clear. The user canexplore the state space by changing the position of the camera in 2D or 3Dspace (effectively changing the Model-View matrix in GPU memory) usingthe mouse and keyboard. The system’s parameters can also be changedduring simulation by dragging sliders on the user interface, which updatesthe location in GPU memory where the corresponding parameter value isstored. This corresponds to glBlendFunc(GL SRC ALPHA, GL ONE) . Appendix B Supplementary Movies • stngpe 2d.mp4 : Simulation of the two-dimensional STN-GPe model,showing the bifurcations that occur as the parameter w ss is increased(refer to main text for details). • stngpe bif.mp4 : Interactive three-dimensional bifurcation diagram forthe STN-GPe model, with the parameter w ss plotted as the third di-mension. • lorenz.mp4 : The Lorenz system visualized in three dimensions. Ini-tially the parameter r has a low value and the system has two stablespirals. After some time, the parameter is changed to its “classical”value of 28, causing the spirals to become unstable and a strange at-tractor to appear. Particle colours as in Figure 3. • hh combined.mp4 : Three Hodgkin-Huxley type neurons coupled by ex-citatory synapses. The position of each particle in three dimensionalspace corresponds to the membrane potential of each neuron. After an23nitial period of unsynchronized firing, the particles all eventually ap-proach of several limit cycles corresponding to different spiking patterns(see main text). 24igure 1: Exploring the two dimensional STN/GPe model in Fireflies un-der variation of w ss . (A)-(D) show screenshots from Fireflies for differentvalues of w ss , with 700,000 particles, half of which are integrated forwardsin time (green) and half of which are integrated backwards (pink). (A) Fullscreenshot of the Fireflies window with w ss = 0. All forward-moving particlesapproach the globally stable spiral (in this still picture it is very difficult tosee the backwards particles, which diverge quickly to infinity). (B) w ss = 4 . w ss = 7 . w ss = 11: An unstable node,saddle (not visible), and stable node. (E) Phase portrait (not generatedby Fireflies) of the STN/GPe system with w ss = 0. All trajectories startfrom the borders of the phase plane and spiral in to the fixed point. Thetrajectories that start in the four corners of the phase plane are shown inpink. 25igure 2: 3D bifurcation diagram of the STN/GPe system under variation ofthe parameter w ss . The simulation contains 8 million particles, divided intoa forwards in time group (green) and a backwards in time group (pink). Thethree visible bifurcations are marked with asterisks, from left to right: fold oflimit cycles, subcritical Andronov-Hopf, saddle node on invariant circle. Notehow the particles on the stable limit cycle become increasingly “bunched”at the top as the cycle approaches the SNIC bifurcation. The asterisks anddirection arrows were added to the image manually and were not generatedby Fireflies. 26igure 3: Exploring the Lorenz equations by varying r . All particles are inte-grated forwards in time, different colours use different initial condition ranges.Asterisks approximately indicate the origin. (A) r = 0 .
5: All particles ap-proach the stable fixed point at the origin. The complex shape is due to theboundaries of the initial conditions ( x ∈ ( − , , y ∈ ( − , , z ∈ (0 , r = 5: The origin is nowa saddle and two new stable fixed points have emerged. Both groups ofparticles have initial conditions near to the origin, but on opposite sidesof its incoming separatrix (green: x, y ∈ (0 , . , z ∈ (0 , . x, y ∈ ( − . , , z ∈ (0 , . r = 15: Trajecto-ries now switch between the two spirals for some time, before settling downand approaching one of them (transient chaos). (C) has initial conditionsnear the saddle, as in (B); these trajectories loop once around one spiral be-fore approaching the other. (D) has initial conditions further away from thesaddle; these trajectories can switch sides repeatedly. (E) r = 25: A strangeattractor. Trajectories starting near the saddle (green and blue) begin byapproaching one of the spirals, before slowly spiralling away and swappingsides chaotically. Another group of particles (bronze) with initial conditionsspanning a large area show the shape of the attractor.27igure 4: A “live” bifurcation diagram for the Lorenz equations. The diagramwas created by setting up a simulation with 8 million particles, each with arandom fixed value for r , projected onto the ( r, y ) plane. The supercriticalpitchfork at r = 1 is clearly visible, as is the appearance of a strange attractorat r ≈ . r = 92. Axes added manually. 28igure 5: A screenshot of Fireflies running the single Hodgkin-Huxley neuronmodel with 500,000 particles. The position of each particle is projected intothe 3D space ( h , m , n ). Particles are coloured according to their positionon 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. 29igure 6: A Hodgkin-Huxley model with three neurons spiking repetitively inresponse to current injection, connected by weak excitatory synapses. (A) Anoverview of the model’s structure. (B) A screenshot of Fireflies running themodel with 500,000 particles. The position of each particle is projected intothe 3D space ( V , V , V3