Pattern phase diagram of spiking neurons on spatial networks
PPattern phase diagram of spiking neurons on spatial networks
Dionysios Georgiadis
Future Resilient Systems, Singapore andDepartment of Management, Technology and Economics,ETH Zurich, Switzerland
Didier Sornette
Department of Management, Technology and Economics,ETH Zurich, Switzerland (Dated: December 24, 2018)We study an abstracted model of neuronal activity via numerical simulation, and report spa-tiotemporal pattern formation and critical like dynamics. A population of pulse coupled, discretised,relaxation oscillators is simulated over networks with varying edge density and spatial embedded-ness. For intermediate edge density and sufficiently strong spatial embeddedness, we observe a novelspatiotemporal pattern in the field of oscillator phases, visually resembling the surface of a frothingliquid. Increasing the edge density results in critical dynamics, with the distribution of neuronalavalanche sizes following a power law with exponent one. Further increase of the edge density resultsin metastable behaviour between pattern formation and synchronisation, before transitioning thesystem entirely into synchrony.
I. INTRODUCTION
Pulse coupled oscillator models (PCO) are definedas populations of relaxation oscillators, interacting ina pulse-like manner over some topology. Such mod-els have provided insight into the phenomenon of spon-taneous synchronisation across fields, be it swarms ofblinking fireflies, pulsing pacemaker heart cells[1], dis-tributed computing systems[2, 3], or traders in financialmarkets [4, 5]. PCOs have been prominent in neurol-ogy, since they can demonstrably capture multiple fea-tures of the rich dynamic behaviour of biological neu-ronal systems, such as: metastability (with the sameconfiguration alternating between synchronous and asyn-chronous behaviour)[6], spatiotemporal pattern forma-tion (in the form of nonlinear waves)[7–11], and criticaldynamics[6, 12, 13].Models used in this literature range from intricate (asfor example the Hodgkin-Huxley model), to relativelysimple (such as Integrate & Fire oscillators). Deville andPeskin [6] studied a population of Discretised Integrate& Fire (DIF) oscillators, arguably the most abstractedmodel of neuronal activity to date. In spite of its simplic-ity, the DIF model with all-to-all stochastic interactionscan successfully capture metastability and synchrony [6],as observed in cortical networks. However, the extendto which the DIF can replicate the phenomenology of itsmore intricate counterparts is not well understood, whichis the motivation for the present paper.Numerical studies have considered the DIF overquenched nontrivial topologies, focusing either on syn-chrony in complex networks[14] or on replicating tem-poral features of biological neuronal networks [13]. No-tably, in [13], populations of cascading spiking neuronsin lab-grown cortical slices were found to follow a power law distribution as well as a specific temporal profile.Both these features were captured by a single simula-tion of a DIF model over a regular lattice, revealing thatcritical-like behaviour is possible in DIF models. Thisobservation has not been studied further - in spite of therecent, sharp interest in models of criticality in neuronalsystems[15–17]. Furthermore, pattern formation in DIFmodels has so far not been studied at all.To address these two issues, we consider the DIF overtopologies ranging between two purposefully picked ex-tremes: a random graph (that has been studied the mostin the literature so far) and a spatial graph (which doesnot disregard the spatial nature of biological systems).Our results showcase that spatial embeddedness endowsthe DIF model with pattern forming behaviour and crit-ical like dynamics, and therefore the unrealistic all-to-allannealing topology drastically restricts the phenomenol-ogy of the DIF. Specifically, for sufficiently low spatialembeddedness, the field defined by the oscillator phasesforms a novel spatiotemporal pattern - visually reminis-cent of the surface of a frothing liquid. Further increas-ing the number of connections between oscillators, wereport metastability between pattern formation and syn-chrony (a behaviour of biological neuronal networks[7]),while the transition itself is characterised by critical-likedynamics. For even higher edge density, near periodicsynchronous firing of all neurons ensues.
II. THE MODEL
Consider N identical oscillators over a biderctional un-weighted graph G . Each node is associated with a binary state a i (1 fired, 0 not fired) and a discrete, non-negative phase φ i . We initialise all oscillator phases uniformly at a r X i v : . [ q - b i o . N C ] D ec random, and all states at zero. We then apply stochasticdrive , by randomly picking a small population of d oscil-lators and increasing their phase by one unit. Variable d is dubbed drive rate . a. Dynamics: Once the phase of an oscillatorreaches the threshold value Θ the oscillator is said to fire :its state is set to one, and the phases of all its neighboursare increased by one. This pulse-like interaction resultsin avalanching events, which we will refer to as cascades .The cascade continues as long as new oscillators fire un-der the influence of their neighbours. Crucially, an os-cillator is only allowed to fire once per cascade - a prop-erty dubbed refractoriness. This is realised in the modelthrough the state a i : an oscillator may only fire if itsstate starts from zero. Thus, once the oscillator fires,and its state is set to one, it is unable to fire again. Thestate remains one until the cascade ceases. Then, all firedoscillators are reset: ( a i , φ i ) ← (0 , , ∀ { i : a i = 1 } . Af-ter resetting all fired oscillators, we resume the stochasticdrive, until a new cascade starts. The number of oscilla-tors that participated in the cascade is the cascade size . b. Topology: The graph of interactions G is a geo-metric random graph[18], furnished with random long-range connections. The random geometric graph is ar-guably the most parsimonious model of spatially embed-ded networks, constituting a canonical choice for the cur-rent study. Given the number of nodes N , the averagedegree E and the fraction of long-to-short range connec-tions R , the topology is assembled in the following way:1. ‘Sprinkle’ N points uniformly at random over a one-by-one square space.2. Find the unconnected pair of points with the min-imum Euclidean distance and connect it. Repeat,until E (1 − R ) edges have been drawn. We willrefer to these edges as short range connections .3. Pick ER random pairs of nodes and connect them,forming the long range connections .For extremely low values of R , we observe a heavily‘meshed’ structure, with a large topological diameter.For R close to one, we retrieve a random graph withzero spatial embeddedness. In the current work, periodicboundary conditions were used in the making of G . c. Additional nomenclature: We measure time indiscrete ‘long’ time, where an additional cascade corre-sponds to an additional ‘long’ time unit. This is to becontrasted with the ‘short’ time, which is associated withthe dynamics within a cascade. A total of n s cascadesare simulated, and the fractional size the t th cascade isdenoted by c t (with c t = 0 if no cascading occured). Ad-ditionally, the phases φ i form a field Φ t = ( φ , . . . , φ N ).Both c t and Φ t vary in time, forming timeseries: C =( c , . . . c n s , ) and Φ = (Φ , . . . , Φ n s ). Both C and Φ aredependent on the parameters E, R, Θ , N . However, forthe sake of notational simplicity, these dependences willremain implicit unless necessary. d. Simulation parameters: Throughout the study,we fix Θ = 5 , d = N/ . Unless otherwise stated,we always simulate 5 · cascades, discarding the first10 events - to ensure that the dynamics have reachedstationarity. This results in n s = 4 · . R, E, N will allbe specified for each experiment.
III. IDENTIFYING REGIMES
For sufficiently high mean degree E , we observe cas-cades of scale O ( N ) over regular time intervals, a prop-erty to which we will refer to as near periodic synchrony .We quantify this behaviour with the help of the metric h ∈ [0 , C ,the time series of cascades on the ‘long’ time process.Values of h near one correspond to asynchrony, while lowvalues of h imply synchrony. For the sake of succinctness,details on this method can be found in Appendix A.For sufficiently low R and over a range of E , we observea spatiotemporal pattern in the Φ field: low phase patches are separated by high-phase ‘fences’ (see first and fourthrows of figure 1). The pattern constantly shifts: cascadesare more likely to occur along the ‘fences’, relaxing the os-cillators and leaving a low phase patch where a fence oncestood. Simultaneously, cascades are unable to propagatethrough large patches, and instead stop in their midst -leaving a ‘fence’ where a patch was before. For the sakeof reference, we dub this spatiotemporal behaviour froth .As depicted in the first and fourth rows of figure 1,the size of the patches increases along with E . As ex-plained above, each patch is the imprint of a past cas-cade - and therefore patch sizes are linked to cascadesizes. Consequently, this behaviour can also be observedin the respective complementary cumulative probabilitydistribution (CCDF) of cascade sizes, presented in thesecond and fourth rows of figure 1. Specifically, as E in-creases, patches are enlarged and the CCDF of cascadesizes extends further towards the right. Eventually, theCCDF forms a power law with exponent one, producingwhat is often referred to as Zipf’s law [19]. At the samepoint, the probability of a global sized patch becomesnonzero. Further increasing E , results in a cascade sizedistribution typical of supercritical dynamics. An exam-ple of supercritical frothing is depicted in the bottomright panel of figure 1. Note the global sized patch, con-nected top-to-bottom, and the supercritical CCDF.To quantify the presence of froth, and the associatedscale of the patches, we consider: s Φ ( (cid:126)λ ) = (cid:104) H t ( (cid:126)λ ) (cid:105) (1)where H t ( (cid:126)λ ) is the spatial Fourier transform of Φ t for thewavelengths (cid:126)λ = [ λ , λ ], and (cid:104) . (cid:105) is the average of multi-ple realisations over time. We exploit the radial symme-FIG. 1: Rows 1 and 4: snapshots of the oscillator phase fields Φ t . The scale of the patterns χ increases along withthe mean degree E . Rows 2 and 5: CCDF of cascade sizes for the simulations of the row above. As χ increases, theCCDF extends towards the left - at first forming a power law and eventually a supercritical distribution (bottomrightmost panel). Rows 3 and 6: blue dots correspond to the temporal average of the spatial spectral power S Φ ( λ )defined in (2). The solid red line is a truncated power law fit of S Φ ( λ ), and the purple (dashed) line is the cornerwavelength χ , as determined by the fitting method described in Appendix B.try of the model by taking the radial mean of s Φ ( (cid:126)λ ): S Φ ( λ ) = 12 π | λ | (cid:90) Ω λ s Φ ( (cid:126)λ ) d Ω λ (2)where Ω λ is a circular shell of radius λ = || (cid:126)λ || .Rows 3 and 6 of figure 1 reveal that frothing is accom-panied by a power law increase of S Φ ( λ ). The increasestarts from a low limit of λ , associated with the averagespatial distance between oscillators neighbouring in thetwo dimensional Euclidean space, and persists up to awavelength where S Φ ( λ ) forms a ‘knee’. The wavelengthassociated with the ‘knee’, dubbed corner wavelength anddenoted by χ , corresponds to the largest spatial scale upto which the froth exhibits a random, self-similar sym-metry, as seen in rows 1 and 4 of figure 1.These observations allow us to numerically quantifythe presence of froth and the size of the patches, by fitting a truncated power law over S Φ ( λ ) (see the red line onrows 3 and 6 of figure 1). The fitting method places thepower law truncation point at a wavelength that servesas an approximation for χ . Additionally, the goodnessof the fit (given by r ∈ [0 , S Φ ( λ )follows a truncated power law. Details on this methodcan be found in Appendix B.The metric h (defined in Appendix A, equation (A2))can be used along with r , to define empirical criteria forthe identification of synchrony and frothing. Concretely,with m h , m r being two threshold constants, we have: h ( E, R, N, Θ) (cid:26) > m h , asynchrony ≤ m h , synchrony (3a) r ( E, R, N, Θ) (cid:26) > m r , frothing up to scale χ ≤ m r , no frothing (3b)FIG. 2: Transitions in the Discretised Integrate & Fire system : h, r and χ √ N (defined in equation (A2), AppendixB §
4, and Section III §
5) as functions of the mean degree E , for a number of system sizes N , with colder colorsindicating larger system sizes (1 . k, . k, k, k, k, k ). Long range connectivity R increases according to ageometric series , from left to right. Low values of h indicate synchrony, and r near one implies the presence ofspatiotemporal patterns. Missing points in the bottom row indicate the absence of pattern formation, for the twolargest R values. Note the transitions i) from high to low values of h and ii) ascending and descending to a plateauof r . The sharpness of both transitions increases with N , demonstrating that larger systems exhibit moreprominent patterns and sharper transitions between macroscopic behaviours. IV. TRANSITION DYNAMICS
The previously proposed criteria (3) require estimatesfor the thresholds m r , m h - which we obtain through theexperiment presented in the current section. R and N are sampled geometrically (as defined in figure 2), while E takes 75 values evenly spaced in [5 , a. Transition to synchrony: In the top row of fig-ure 2 we observe a sharp transition, from high to low h , revealing that the system suddenly moves from asyn-chrony to synchrony beyond a value of E . As R in-creases, the transition shifts to smaller E values, reveal-ing that spatial embeddedness delays the onset of syn-chrony. The sharpness of the transition increases alongwith system size, demonstrating that this transition willstill be present - and even more prominent - in the ther-modynamic limit. b. Frothing delays synchrony: The second row in fig-ure 2 depicts a clear plateau of the r metric, implying thepresence of frothing dynamics over a range of E values.The frothing regime is interposed between asynchronyand synchrony, with its width decreasing as R increases.Therefore, while high R systems enter synchrony, low R systems froth instead - implying that frothing is themechanism that delays synchrony. The presence of thefrothing regime does not depend on size, since the plateauof r persists - and even widens - along with N . c. Scaling of χ : The bottom row of figure 2 showsthat the corner frequency χ increases with N . Specifi-cally, for the range of values in this study (N from 1 . k to 40 k ) plotting χ √ N makes simulation results for all N to collapse into a single universal curve, revealing thepresence of a scaling law. Finally, for all panels in thebottom row of figure 2, the peak of χ coincides with theend of the r plateau, verifying that frothing patternsbecome the most prominent on the verge of synchrony. V. EMPIRICAL REGIME DIAGRAM
To numerically investigate the macroscopic behaviourof the model over the
E, R space, we fix N = 10 , andvary R over 70 evenly spaced values in [6 , E takes 30geometrically-spaced values in [10 − , m h , m r . Visual inspectionof the results in figure 2 reveals that m h = 5 10 − and m r = 0 . • Regime I:
For low connection density the systemexhibits local cascades, with no frothing. Theliterature[6, 20] refers to this regime as asynchrony . • Regime II:
For low long range connectivity, andover a mid-range of E , frothing appears in the Φfield. We refer to this regime as frothing regime . • Regime III:
Starting from froth and sufficiently in-creasing E results in a regime where we intermit-tently observe the phenomenology of regimes II andIV. The CCDF of cascade sizes is characteristic ofa supercritical system, with slower than power lawdecay, and global-sized cascades appearing regu-larly. This regime is dubbed metastable regime . • Regime IV:
For high connection density, the sys-tem undergoes a discontinuous limit cycle. As theaverage phase increases with time, cascades remainlocal in scale, until a global sized cascade occurs -resulting in the relaxation of the phase field. Then,the buildup of phase synchronisation begins anew.We refer to this regime as synchrony .FIG. 3: Empirical regime diagram of the DiscretisedIntegrate & Fire model over the space of degree densityand long range connectivity (
E, R ). We observe fourregimes: asynchrony (I), pattern formation (II),synchrony (IV), and metastability between patternformation and synchrony (III). The regime boundariesare drawn based on the criteria defined in (3) . VI. DISCUSSION
We have numerically investigated the behaviour of dis-cretised Integrate & Fire oscillators, with slow stochasticdrive, over spatial graphs. Remarkably, when placed overspatial topologies these models give rise to novel, nontriv-ial dynamics: spatiotemporal patterns form, where largeclusters of relaxed nodes act as natural barriers againstcascading, while jagged strips of near-firing nodes facili-tate long distance propagation of cascades. This patterncan give rise to critical dynamics with cascading eventssizes following a power law with exponent one. Furtherincreasing the number of edges results in metastabilitybetween pattern formation and synchrony, and eventu-ally drives the model into synchrony.Our findings show that taking into considerationthe fundamentally spatial character of neuronal net-works drastically extends the phenomenological overlap between discretised Integrate & Fire and more intri-cate neuronal models. This observation provides solidgrounds for the usage of the discretised Integrate & Firemodel in the study of pattern formation and critical dy-namics. Further numerical studies could shed light on therole of additional neuronal properties (leakiness, nonlin-ear response curves, and inhibition) on pattern formationand criticality in neuronal systems.In spite of the extensive simulations in the currentstudy, the exact mathematical nature of the regime tran-sitions of the model is not well understood - and couldconstitute the subject of future works. Additionally, tothe best of the authors’ knowledge, frothing has not beenso far empirically observed in physical systems. This issurprising considering the generality of the conditions un-der which frothing arises. A possible explanation couldbe that detecting froth requires knowledge of the phase ,an attribute of oscillators that is typically less prominentand harder to measure than their state . In any case, fur-ther works are warranted, focusing on the empirical de-tection of frothing in Pulse Coupled Oscillator systems.One of the authors (DG) would like to thank S.C. Leraand S.O.P. Blume for stimulating discussions, and is sup-ported by the Future Resilient Systems at the Singapore-ETH Centre (SEC), a collaboration between ETH Zurichand National Research Foundation (NRF) of Singapore(FI 370074011) - under the auspices of the NRF’s Cam-pus for Research Excellence and Technological Enterprise(CREATE) programme.
Appendix A: Quantifying synchrony
The near-periodic behaviour of the time series of cas-cade sizes C can be detected in the frequency domain,where peaks appear in the power spectrum. The con-centration of the spectral power around these peaks isquantified using the Herfindahl-Hirschman index. Sincethe data are discrete, we will be using the discrete Fouriertransform. Let P C ( f ) be the temporal spectral power of C for frequency f :ˆ h = (cid:88) f (cid:18) P C ( f ) N (cid:19) , N = (cid:88) f P C ( f ) (A1) h = ˆ h − n s − n s (A2)The ˆ h metric takes values in [0 ,
1] and quantifies the con-centration of P C ( f ) around a few frequencies. Duringasynchrony, the only significantly contributing frequencyis the 0th, resulting in near zero values of P C ( f ) for f > h near one. In contrast, during synchrony,higher harmonics carry considerable amount of power re-sulting in a lower h index. Discerning between synchronyand asynchrony does require a threshold value of h , whichwe determine empirically in this study to be equal to5 · − (see Section V § Appendix B: Quantifying spatiotemporal patternformation
We observe that in the case of frothing patterns, S Φ ( λ )follows a power law increase, over a band of wavelengths.The upper limit of the band χ is associated with thelargest cell size of the frothing pattern. The lower limitof the band is proportional to the size of the mesh used toestimate the Φ field. In the current study, √ N oscillatorsare placed along each dimension of the two-dimensionalEuclidean space, resulting in an lower wavelength limitof 4 π/ √ N .Higher frequency components also reside in the Φfield, for example due to the randomness in the place-ment of the oscillators. These higher frequency compo-nents may introduce noise in the spatial spectrum of Φ,through a processes known as aliasing . Specifically inthe case of data produced by processes with power lawspectra, aliasing results in the measured spectrum pro-gressively resembling white noise as we move to smallerwavelengths[21]. As a treatment, we ignore the values ofthe measured spectrum for small wavelengths, by dou-bling the lower limit derived in the previous paragraphto λ min = 8 π/ √ N .To determine the corner wavelength χ , we consider the frequency response function of a linear low pass filter: g ( λ ; p , p , p , p ) = p (cid:113) (cid:0) λ/p (cid:1) − p + p (B1)and fit it to S Φ ( λ ), according to the following equation:( p ∗ , p ∗ , p ∗ , p ∗ ) =argmin p ,p ,p ,p (cid:88) λ ≥ λ min (cid:18) S Φ ( λ ) − g ( λ ; p , p , p , p ) S Φ ( λ ) (cid:19) (B2)where p , p , p , p are fitting parameters. Equation (B1)describes a power law increase with exponent p , upto the wavelength p where the function forms a visual‘knee’. From that point onwards, the value of (B1) re-mains nearly constant at p + p . For the sake of illustra-tion, and to enable comparison between different S Φ ( λ )curves, we position the truncation point χ at p ∗ .Initial solutions to the problem in (B2) were ob-tained via the particle swarm method, and refined witha Levenberg-Marquardt local search. Examples of thefitted g ( λ ; p ∗ , p ∗ , p ∗ , p ∗ ) are depicted in rows 3 and 6 ofFigure 1 (red solid line), along with the corresponding S Φ ( λ ) values. The quality of the fits can be assessedvia the r-squared metric, to which we will refer to as r .Specifically, r is used to determine whether a simulationexhibits frothing: high value of r implies that equationB1 provides a good fit, providing evidence in support ofthe presence of froth in the Φ field. [1] R. E. Mirollo and S. H. Strogatz, SIAM Journal on Ap-plied Mathematics , 1645 (1990).[2] A. V. Proskurnikov and M. Cao, IEEE Transactions onAutomatic Control , 5873 (2017).[3] K. Konishi and H. Kokame, Chaos: An InterdisciplinaryJournal of Nonlinear Science , 033132 (2008).[4] C. M. Wray and S. R. Bishop, Scientific Reports , 6355(2014).[5] C. M. Wray and S. R. Bishop, PloS One , e0151790(2016).[6] R. L. DeVille and C. S. Peskin, Bulletin of MathematicalBiology , 1608 (2008).[7] X. Huang, W. C. Troy, Q. Yang, H. Ma, C. R. Laing,S. J. Schiff, and J.-Y. Wu, Journal of Neuroscience ,9897 (2004).[8] H. R. Wilson and J. D. Cowan, Biophysical Journal ,1 (1972).[9] R. L. DeVille, E. Vanden-Eijnden, and C. B. Muratov,Physical Review E , 031105 (2005).[10] G. B. Ermentrout and D. Kleinfeld, Neuron , 33(2001).[11] X. Guardiola and A. D´ıaz-Guilera, Phys. Rev. E , 3626 (1999).[12] S. Bottani, Physical Review Letters , 4189 (1995).[13] N. Friedman, S. Ito, B. A. Brinkman, M. Shimono, R. L.DeVille, K. A. Dahmen, J. M. Beggs, and T. C. Butler,Physical Review Letters , 208102 (2012).[14] R. L. DeVille and C. S. Peskin, Bulletin of MathematicalBiology , 769 (2012).[15] J. M. Beggs and D. Plenz, Journal of Neuroscience ,11167 (2003).[16] O. Kinouchi and M. Copelli, Nature Physics , 348(2006).[17] J. M. Beggs, Philosophical Transactions of the Royal So-ciety of London A: Mathematical, Physical and Engineer-ing Sciences , 329 (2008).[18] M. Penrose et al. , Random geometric graphs , 5 (Oxforduniversity press, 2003).[19] A. Saichev, Y. Malevergne, and D. Sornette,
Theory ofZipf’s Law and Beyond (Springer, 2010).[20] R. L. DeVille, C. S. Peskin, and J. H. Spencer, Mathe-matical Modelling of Natural Phenomena , 26 (2010).[21] J. W. Kirchner, Physical Review E71