Pattern phase diagram of spiking neurons on spatial networks

TL;DR

This paper investigates how spiking neurons on spatial networks form complex patterns and exhibit critical dynamics, revealing a transition from pattern formation to synchrony as network connectivity increases.

Contribution

It introduces a novel spatiotemporal pattern in neuronal phase fields and characterizes the transition to critical and synchronized states based on network parameters.

Findings

Discovery of a frothing liquid-like spatiotemporal pattern.

Neuronal avalanche sizes follow a power law with exponent one.

Identification of metastable behavior before full synchronization.

Abstract

We study an abstracted model of neuronal activity via numerical simulation, and report spatiotemporal 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 novel spatiotemporal pattern in the field of oscillator phases, visually resembling the surface of a frothing liquid. Increasing the edge density results in critical dynamics, with the distribution of neuronal avalanche sizes following a power law with exponent one. Further increase of the edge density results in metastable behaviour between pattern formation and synchronisation, before transitioning the system entirely into synchrony.

Equations

$$s_{\Phi }(\vec{\lambda })=⟨H_t^2(\vec{\lambda })⟩$$

$${\mathcal{S}}_{\Phi }(\lambda )=\frac{1}{2\pi |\lambda |}{\int }_{{\Omega }_{\lambda }}{s}_{\Phi }(\vec{\lambda })d{\Omega }_{\lambda }$$

$$h(E,R,N,\Theta )$$

$$r^2(E,R,N,\Theta )$$

$$\hat{h}=\sum _f{\left(\frac{{\mathcal{P}}_C(f)}{\mathcal{N}}\right)}^2,\mathcal{N}=\sum _f{\mathcal{P}}_C(f)$$

$$h=\frac{\hat{h}-n_s}{1-n_s}$$

$$g(\lambda;p_{1},p_{2},p_{3},p_{4})=\frac{p_{1}}{\sqrt{1+\big{(}\lambda/p_{3}\big{)}^{-2p_{4}}}}+p_{2}$$

$$\begin{gathered} (p_1^{\ast },p_2^{\ast },p_3^{\ast },p_4^{\ast })= \\ \underset{\begin{array}{c}{p}_1,p_2,\\ p_3,p_4\end{array}}{\mathrm{argmin}}\sum _{\lambda \ge {\lambda }_{\text{min}}}{\left(\frac{{\mathcal{S}}_{\Phi }(\lambda )-g(\lambda ;p_1,p_2,p_3,p_4)}{{\mathcal{S}}_{\Phi }(\lambda )}\right)}^2 \end{gathered}$$

Full text

I Introduction

Pulse coupled oscillator models (PCO) are defined as populations of relaxation oscillators, interacting in a pulse-like manner over some topology. Such models have provided insight into the phenomenon of spontaneous synchronisation across fields, be it swarms of blinking fireflies, pulsing pacemaker heart cellsMirollo and Strogatz (1990), distributed computing systemsProskurnikov and Cao (2017); Konishi and Kokame (2008), or traders in financial markets Wray and Bishop (2014, 2016). PCOs have been prominent in neurology, since they can demonstrably capture multiple features of the rich dynamic behaviour of biological neuronal systems, such as: metastability (with the same configuration alternating between synchronous and asynchronous behaviour)DeVille and Peskin (2008), spatiotemporal pattern formation (in the form of nonlinear waves)Huang et al. (2004); Wilson and Cowan (1972); DeVille et al. (2005); Ermentrout and Kleinfeld (2001); Guardiola and Díaz-Guilera (1999), and critical dynamicsDeVille and Peskin (2008); Bottani (1995); Friedman et al. (2012).

Models used in this literature range from intricate (as for example the Hodgkin-Huxley model), to relatively simple (such as Integrate & Fire oscillators). Deville and Peskin DeVille and Peskin (2008) studied a population of Discretised Integrate & Fire (DIF) oscillators, arguably the most abstracted model of neuronal activity to date. In spite of its simplicity, the DIF model with all-to-all stochastic interactions can successfully capture metastability and synchrony DeVille and Peskin (2008), as observed in cortical networks. However, the extend to which the DIF can replicate the phenomenology of its more intricate counterparts is not well understood, which is the motivation for the present paper.

Numerical studies have considered the DIF over quenched nontrivial topologies, focusing either on synchrony in complex networksDeVille and Peskin (2012) or on replicating temporal features of biological neuronal networks Friedman et al. (2012). Notably, in Friedman et al. (2012), populations of cascading spiking neurons in 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 simulation of a DIF model over a regular lattice, revealing that critical-like behaviour is possible in DIF models. This observation has not been studied further - in spite of the recent, sharp interest in models of criticality in neuronal systemsBeggs and Plenz (2003); Kinouchi and Copelli (2006); Beggs (2008). Furthermore, pattern formation in DIF models has so far not been studied at all.

To address these two issues, we consider the DIF over topologies ranging between two purposefully picked extremes: a random graph (that has been studied the most in the literature so far) and a spatial graph (which does not disregard the spatial nature of biological systems). Our results showcase that spatial embeddedness endows the DIF model with pattern forming behaviour and critical like dynamics, and therefore the unrealistic all-to-all annealing topology drastically restricts the phenomenology of the DIF. Specifically, for sufficiently low spatial embeddedness, the field defined by the oscillator phases forms a novel spatiotemporal pattern - visually reminiscent of the surface of a frothing liquid. Further increasing the number of connections between oscillators, we report metastability between pattern formation and synchrony (a behaviour of biological neuronal networksHuang et al. (2004)), while the transition itself is characterised by critical-like dynamics. For even higher edge density, near periodic synchronous firing of all neurons ensues.

II The model

Consider $N$ identical oscillators over a biderctional unweighted graph $\mathcal{G}$. Each node is associated with a binary state $a_{i}$ ($1$ fired, [math] not fired) and a discrete, non-negative phase $\phi_{i}$. We initialise all oscillator phases uniformly at random, and all states at zero. We then apply stochastic drive, by randomly picking a small population of $d$ oscillators and increasing their phase by one unit. Variable $d$ is dubbed drive rate.

Dynamics:

Once the phase of an oscillator reaches the threshold value $\Theta$ the oscillator is said to fire: its state is set to one, and the phases of all its neighbours are increased by one. This pulse-like interaction results in avalanching events, which we will refer to as cascades. The cascade continues as long as new oscillators fire under the influence of their neighbours. Crucially, an oscillator is only allowed to fire once per cascade - a property dubbed refractoriness. This is realised in the model through the state $a_{i}$: an oscillator may only fire if its state starts from zero. Thus, once the oscillator fires, and its state is set to one, it is unable to fire again. The state remains one until the cascade ceases. Then, all fired oscillators are reset: $(a_{i},\phi_{i})\leftarrow(0,0),~{}\forall~{}\{i:a_{i}=1\}$. After resetting all fired oscillators, we resume the stochastic drive, until a new cascade starts. The number of oscillators that participated in the cascade is the cascade size.

Topology:

The graph of interactions $\mathcal{G}$ is a geometric random graphPenrose et al. (2003), furnished with random long-range connections. The random geometric graph is arguably the most parsimonious model of spatially embedded networks, constituting a canonical choice for the current study. Given the number of nodes $N$, the average degree $E$ and the fraction of long-to-short range connections $R$, the topology is assembled in the following way:

‘Sprinkle’ $N$ points uniformly at random over a one-by-one square space. 2. 2.

Find the unconnected pair of points with the minimum Euclidean distance and connect it. Repeat, until $E(1-R)$ edges have been drawn. We will refer to these edges as short range connections. 3. 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 with zero spatial embeddedness. In the current work, periodic boundary conditions were used in the making of $\mathcal{G}$.

Additional nomenclature:

We measure time in discrete ‘long’ time, where an additional cascade corresponds to an additional ‘long’ time unit. This is to be contrasted with the ‘short’ time, which is associated with the dynamics within a cascade. A total of $n_{s}$ cascades are simulated, and the fractional size the $t$th cascade is denoted by $c_{t}$ (with $c_{t}=0$ if no cascading occured). Additionally, the phases $\phi_{i}$ form a field $\Phi_{t}=(\phi_{0},\dots,\phi_{N})$. Both $c_{t}$ and $\Phi_{t}$ vary in time, forming timeseries: $C=(c_{0},\dots c_{n_{s}},)$ and $\Phi=(\Phi_{0},\dots,\Phi_{n_{s}})$. Both $C$ and $\Phi$ are dependent on the parameters $E,R,\Theta,N$. However, for the sake of notational simplicity, these dependences will remain implicit unless necessary.

Simulation parameters:

Throughout the study, we fix $\Theta=5,~{}d=N/{10^{3}}$. Unless otherwise stated, we always simulate $5\cdot 10^{4}$ cascades, discarding the first $10^{4}$ events - to ensure that the dynamics have reached stationarity. This results in $n_{s}=4\cdot 10^{4}$. $R,E,N$ will all be specified for each experiment.

III Identifying regimes

For sufficiently high mean degree $E$, we observe cascades of scale $\mathcal{O}(N)$ over regular time intervals, a property to which we will refer to as near periodic synchrony. We quantify this behaviour with the help of the metric $h\in[0,1]$, which is formally defined as the normalised Herfindahl index of the temporal power spectrum of $C$, the time series of cascades on the ‘long’ time process. Values of $h$ near one correspond to asynchrony, while low values 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 observe a spatiotemporal pattern in the $\Phi$ field: low phase patches are separated by high-phase ‘fences’ (see first and fourth rows of figure 1). The pattern constantly shifts: cascades are more likely to occur along the ‘fences’, relaxing the oscillators and leaving a low phase patch where a fence once stood. Simultaneously, cascades are unable to propagate through large patches, and instead stop in their midst - leaving a ‘fence’ where a patch was before. For the sake of 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 explained above, each patch is the imprint of a past cascade - and therefore patch sizes are linked to cascade sizes. Consequently, this behaviour can also be observed in the respective complementary cumulative probability distribution (CCDF) of cascade sizes, presented in the second and fourth rows of figure 1. Specifically, as $E$ increases, patches are enlarged and the CCDF of cascade sizes extends further towards the right. Eventually, the CCDF forms a power law with exponent one, producing what is often referred to as Zipf’s law Saichev et al. (2010). At the same point, the probability of a global sized patch becomes nonzero. Further increasing $E$, results in a cascade size distribution typical of supercritical dynamics. An example of supercritical frothing is depicted in the bottom right panel of figure 1. Note the global sized patch, connected top-to-bottom, and the supercritical CCDF.

To quantify the presence of froth, and the associated scale of the patches, we consider:

[TABLE]

where $H_{t}(\vec{\lambda})$ is the spatial Fourier transform of $\Phi_{t}$ for the wavelengths $\vec{\lambda}=[\lambda_{1},\lambda_{2}]$, and $\langle.\rangle$ is the average of multiple realisations over time. We exploit the radial symmetry of the model by taking the radial mean of $s_{\Phi}(\vec{\lambda})$:

[TABLE]

where $\Omega_{\lambda}$ is a circular shell of radius $\lambda=||\vec{\lambda}||$.

Rows 3 and 6 of figure 1 reveal that frothing is accompanied by a power law increase of $\mathcal{S}_{\Phi}(\lambda)$. The increase starts from a low limit of $\lambda$, associated with the average spatial distance between oscillators neighbouring in the two dimensional Euclidean space, and persists up to a wavelength where $\mathcal{S}_{\Phi}(\lambda)$ forms a ‘knee’. The wavelength associated with the ‘knee’, dubbed corner wavelength and denoted by $\chi$, corresponds to the largest spatial scale up to which the froth exhibits a random, self-similar symmetry, as seen in rows 1 and 4 of figure 1.

These observations allow us to numerically quantify the presence of froth and the size of the patches, by fitting a truncated power law over $\mathcal{S}_{\Phi}(\lambda)$ (see the red line on rows 3 and 6 of figure 1). The fitting method places the power law truncation point at a wavelength that serves as an approximation for $\chi$. Additionally, the goodness of the fit (given by $r^{2}\in[0,1]$) quantifies how well $\mathcal{S}_{\Phi}(\lambda)$ follows a truncated power law. Details on this method can be found in Appendix B.

The metric $h$ (defined in Appendix A, equation (5)) can be used along with $r^{2}$, to define empirical criteria for the identification of synchrony and frothing. Concretely, with $m_{h},m_{r^{2}}$ being two threshold constants, we have:

[TABLE]

IV Transition dynamics

The previously proposed criteria (3) require estimates for the thresholds $m_{r^{2}},m_{h}$ - which we obtain through the experiment 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,20]$. The results, shown in figure 2, reveal abrupt transitions in the macroscopic behaviour of the model, allowing us to draw three conclusions:

Transition to synchrony:

In the top row of figure 2 we observe a sharp transition, from high to low $h$, revealing that the system suddenly moves from asynchrony to synchrony beyond a value of $E$. As $R$ increases, the transition shifts to smaller $E$ values, revealing that spatial embeddedness delays the onset of synchrony. The sharpness of the transition increases along with system size, demonstrating that this transition will still be present - and even more prominent - in the thermodynamic limit.

Frothing delays synchrony:

The second row in figure 2 depicts a clear plateau of the $r^{2}$ metric, implying the presence of frothing dynamics over a range of $E$ values. The frothing regime is interposed between asynchrony and synchrony, with its width decreasing as $R$ increases. Therefore, while high $R$ systems enter synchrony, low $R$ systems froth instead - implying that frothing is the mechanism that delays synchrony. The presence of the frothing regime does not depend on size, since the plateau of $r^{2}$ persists - and even widens - along with $N$.

Scaling of $\chi$:

The bottom row of figure 2 shows that the corner frequency $\chi$ increases with $N$. Specifically, for the range of values in this study (N from $1.25k$ to $40k$) plotting $\chi\sqrt{N}$ makes simulation results for all $N$ to collapse into a single universal curve, revealing the presence of a scaling law. Finally, for all panels in the bottom row of figure 2, the peak of $\chi$ coincides with the end of the $r^{2}$ plateau, verifying that frothing patterns become the most prominent on the verge of synchrony.

V Empirical regime diagram

To numerically investigate the macroscopic behaviour of the model over the $E,R$ space, we fix $N=10^{4}$, and vary $R$ over $70$ evenly spaced values in $[6,20]$. $E$ takes $30$ geometrically-spaced values in $[10^{-3},1]$. In order to use the criteria (3), we need to set $m_{h},m_{r^{2}}$. Visual inspection of the results in figure 2 reveals that $m_{h}=5~{}10^{-2}$ and $m_{r^{2}}=0.9$ allow criteria (3) to separate the macroscopic regimes adequately well for illustrative purposes. The resulting empirical regime diagram is depicted in figure 3, revealing four regimes:

• •

Regime I: For low connection density the system exhibits local cascades, with no frothing. The literatureDeVille and Peskin (2008); DeVille et al. (2010) refers to this regime as asynchrony.

• •

Regime II: For low long range connectivity, and over a mid-range of $E$, frothing appears in the $\Phi$ field. We refer to this regime as frothing regime.

• •

Regime III: Starting from froth and sufficiently increasing $E$ results in a regime where we intermittently observe the phenomenology of regimes II and IV. The CCDF of cascade sizes is characteristic of a supercritical system, with slower than power law decay, and global-sized cascades appearing regularly. This regime is dubbed metastable regime.

• •

Regime IV: For high connection density, the system undergoes a discontinuous limit cycle. As the average phase increases with time, cascades remain local 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.

VI Discussion

We have numerically investigated the behaviour of discretised Integrate & Fire oscillators, with slow stochastic drive, over spatial graphs. Remarkably, when placed over spatial topologies these models give rise to novel, nontrivial dynamics: spatiotemporal patterns form, where large clusters of relaxed nodes act as natural barriers against cascading, while jagged strips of near-firing nodes facilitate long distance propagation of cascades. This pattern can give rise to critical dynamics with cascading events sizes following a power law with exponent one. Further increasing the number of edges results in metastability between pattern formation and synchrony, and eventually drives the model into synchrony.

Our findings show that taking into consideration the fundamentally spatial character of neuronal networks drastically extends the phenomenological overlap between discretised Integrate & Fire and more intricate neuronal models. This observation provides solid grounds for the usage of the discretised Integrate & Fire model in the study of pattern formation and critical dynamics. Further numerical studies could shed light on the role of additional neuronal properties (leakiness, nonlinear response curves, and inhibition) on pattern formation and criticality in neuronal systems.

In spite of the extensive simulations in the current study, the exact mathematical nature of the regime transitions of the model is not well understood - and could constitute the subject of future works. Additionally, to the best of the authors’ knowledge, frothing has not been so far empirically observed in physical systems. This is surprising considering the generality of the conditions under which frothing arises. A possible explanation could be that detecting froth requires knowledge of the phase, an attribute of oscillators that is typically less prominent and harder to measure than their state. In any case, further works are warranted, focusing on the empirical detection of frothing in Pulse Coupled Oscillator systems.

One of the authors (DG) would like to thank S.C. Lera and S.O.P. Blume for stimulating discussions, and is supported by the Future Resilient Systems at the Singapore-ETH Centre (SEC), a collaboration between ETH Zurich and National Research Foundation (NRF) of Singapore (FI 370074011) - under the auspices of the NRF’s Campus for Research Excellence and Technological Enterprise (CREATE) programme.

Appendix A Quantifying synchrony

The near-periodic behaviour of the time series of cascade sizes $C$ can be detected in the frequency domain, where peaks appear in the power spectrum. The concentration of the spectral power around these peaks is quantified using the Herfindahl-Hirschman index. Since the data are discrete, we will be using the discrete Fourier transform. Let $\mathcal{P}_{C}(f)$ be the temporal spectral power of $C$ for frequency $f$:

[TABLE]

The $\hat{h}$ metric takes values in $[0,1]$ and quantifies the concentration of $\mathcal{P}_{C}(f)$ around a few frequencies. During asynchrony, the only significantly contributing frequency is the [math]th, resulting in near zero values of $\mathcal{P}_{C}(f)$ for $f>0$ and therefore $h$ near one. In contrast, during synchrony, higher harmonics carry considerable amount of power resulting in a lower $h$ index. Discerning between synchrony and asynchrony does require a threshold value of $h$, which we determine empirically in this study to be equal to $5\cdot 10^{-2}$ (see Section V §1 for details).

Appendix B Quantifying spatiotemporal pattern formation

We observe that in the case of frothing patterns, $\mathcal{S}_{\Phi}(\lambda)$ follows a power law increase, over a band of wavelengths. The upper limit of the band $\chi$ is associated with the largest cell size of the frothing pattern. The lower limit of the band is proportional to the size of the mesh used to estimate the $\Phi$ field. In the current study, $\sqrt{N}$ oscillators are placed along each dimension of the two-dimensional Euclidean space, resulting in an lower wavelength limit of $4\pi/\sqrt{N}$.

Higher frequency components also reside in the $\Phi$ field, for example due to the randomness in the placement of the oscillators. These higher frequency components may introduce noise in the spatial spectrum of $\Phi$, through a processes known as aliasing. Specifically in the case of data produced by processes with power law spectra, aliasing results in the measured spectrum progressively resembling white noise as we move to smaller wavelengthsKirchner (2005). As a treatment, we ignore the values of the measured spectrum for small wavelengths, by doubling the lower limit derived in the previous paragraph to $\lambda_{\text{min}}=8\pi/\sqrt{N}$.

To determine the corner wavelength $\chi$, we consider the frequency response function of a linear low pass filter:

[TABLE]

and fit it to $\mathcal{S}_{\Phi}(\lambda)$, according to the following equation:

[TABLE]

where $p_{1},p_{2},p_{3},p_{4}$ are fitting parameters. Equation (6) describes a power law increase with exponent $p_{4}$, up to the wavelength $p_{3}$ where the function forms a visual ‘knee’. From that point onwards, the value of (6) remains nearly constant at $p_{1}+p_{2}$. For the sake of illustration, and to enable comparison between different $\mathcal{S}_{\Phi}(\lambda)$ curves, we position the truncation point $\chi$ at $p_{3}^{*}$.

Initial solutions to the problem in (7) were obtained via the particle swarm method, and refined with a Levenberg-Marquardt local search. Examples of the fitted $g(\lambda;p_{1}^{*},p_{2}^{*},p_{3}^{*},p_{4}^{*})$ are depicted in rows 3 and 6 of Figure 1 (red solid line), along with the corresponding $\mathcal{S}_{\Phi}(\lambda)$ values. The quality of the fits can be assessed via the r-squared metric, to which we will refer to as $r^{2}$. Specifically, $r^{2}$ is used to determine whether a simulation exhibits frothing: high value of $r^{2}$ implies that equation 6 provides a good fit, providing evidence in support of the presence of froth in the $\Phi$ field.