Spiral wave dynamics in a neuronal network model

Spiral waves are spatial-temporal patterns that can emerge in different systems as heart tissues, chemical oscillators, ecological networks and the brain. These waves have been identified in the neocortex of turtles, rats, and humans, particularly during sleep-like states. Although their functions in cognitive activities remain until now poorly understood, these patterns are related to cortical activity modulation and contribute to cortical processing. In this work, we construct a neuronal network layer based on the spatial distribution of pyramidal neurons. Our main goal is to investigate how local connectivity and coupling strength are associated with the emergence of spiral waves. Therefore, we propose a trustworthy method capable of detecting different wave patterns, based on local and global phase order parameters. As a result, we find that the range of connection radius (R) plays a crucial role in the appearance of spiral waves. For R<20 {\mu}m, only asynchronous activity is observed due to small number of connections. The coupling strength (gsyn ) greatly influences the pattern transitions for higher R, where spikes and bursts firing patterns can be observed in spiral and non-spiral waves. Finally, we show that for some values of R and gsyn bistable states of wave patterns are obtained.


Introduction
It has been identified astonishing spatiotemporal patterns which can be seen in 2D and 3D networks with non-local interactions, named spiral wave chimera states.It consists of synchronous oscillators which rotate around a desynchronized core, which is called phase singularity (PS) [1].This pattern is seen in a variety of fields.In heart tissues [2,3,4] spiral waves have been extensively studied due their correlation with some heart issues, such as ventricular tachycardia and ventricular fibrillation.Davidenko et al. [5] demonstrated the emergence of spiral waves in sheep and dog epicardial muscle by means of potentiometric dye with charge-coupled device imaging technology.The authors showed that sometimes the phase singularity drifted away from its original position and dissipated in the tissue border, this drift was associated with a Doppler shift.Supression of spiral waves in the atrial cardiomyocytes can be done by means of optogenetics [6].Spiral waves are able to emerge in ecological network composed of diffusible prey-predator species locally coupled [7].Totz et al. [8] observed the emergence of spiral wave chimeras in large populations of coupled chemical oscillators, they have studied the motion and splitting of the asynchronous core.
The chimera is defined as the coexistence of coherent and incoherent domains when similar elements interact [9].The notion of chimera states has been discussed in the last decade [10].Majhi et al. [11] reviewed important aspects of such dynamics and highlighted the different types, as well as their relevance in a biophysical context.Chimera states play an important role in empirical brain networks, specially in the cerebral cortex [12].It has been found in networks composed of coupled neurons, such as integrate-and-fire [14,13], Hodgkin-Huxley [9], Hindmarsh-Rose [15], FitzHugh-Nagumo [17,16], Morris-Lecar [18], Nekorkin maps [19] and others.The spatiotemporal patterns, in which coherence coexists with incoherence, have been reported in local [20,21] and non-local network configurations [22,23].
Many studies have improved the understanding of chimera states in neuronal networks.The parameter regions for chimera states in a two coupling configuration of identical thermally sensitive Hodgkin-Huxley neurons was characterized in Ref. [9].Majhi [24] investigated chimera states in uncoupled neurons induced by multilayer interactions.They found that the competition between electrical and chemical synapses can lead to chimeras.Hövel et al. [25] studied the coexistence of chimera and traveling waves, as well as multi-stability, in the parameter space of a neuronal system of excitability type I.Chemical and electrical synapses can affect the features of chimera states and traveling waves [26].The coexistence of patterns can emerge due to electrical and chemical coupling in a multi-weighted of the memristive Fitzhugh-Nagumo network [16].It was demonstrated the importance of such interactions to obtain synchronous and desynchronous behaviours coexisting [27].Chimera states were also reported in time-delay networks [28,29].
Neuronal networks can exhibit spiral wave spatial patterns when in a sleep-like state [31,30].It is able to modify ongoing activity in the cortex [30,32].Local excitatory interactions are of great importance in the generation of spiral waves [30].The role of spiral waves in the neocortex is not yet well understood, thus simulations and experimental methods are crucial to better understand spiral wave dynamics.The emergence of spiral waves in computational simulations has been observed in different models throughout the time.Ma et al. [33] studied the robustness of the spiral wave when noise is applied to ion channels of the standard Hodking-Huxley model.They showed that the spiral wave does not sustain its activity, whether the noise increases or the fraction of active channels is below a threshold.The transmission delay enhances the coherence of spiral waves in a noisy network [34].Spiral waves can be seen in non-locally coupled maps for some coupling parameters [35].Santos et.al. [36] demonstrated that chimera spiral waves with multi-cores appear in regular and fractal networks.Their study showed that the initial condition of the adaptation current has a crucial role in the emergence of spiral waves.
The formation of spiral waves are not only observed in neuronal simulations, but also in experimental data [37,39,38].The effects of spiral waves in cognitive activity are not well defined, although it influences the frequencies [30].Prechtl et al. [38] used visual stimuli to register the appearance of spiral waves in the turtle cortex.Furthermore, Huang et al. [31] captured the emergence of spiral waves in rat cortex and its influence in coordinating oscillations of neuronal populations.Spiral waves have been identified in the neocortex during sleep-like states and pharmacologically induced oscillations [30].It is theorized that such spatial pattern plays an important role in organizing and modulating healthy and pathological brain activities, for instance, epilepsy [31,30,38].
In this work, we propose a neuronal network composed of pyramidal cells, that are identified as excitatory neurons.Pyramidal neurons are the most common cells in the brain, and the name is due to its characteristic shape [40].This neuron type is also the main one of the excitatory neuron family in the brain [40].Our network model is inspired by a two dimensional slice of pyramidal layer of the hippocampus [41].The model considered to describe a single neuron membrane potential is the adaptive exponential integrate-and-fire (aEIF) [42].Besides the simplicity of this model and its low computational cost, it exhibits a great biophysical accuracy [43,44].Due to the lack of reliable tools able to characterize chimera spiral waves, we propose a method to identify chimeras spiral wave, and analyze their dependence on the coupling radius (R) and the synapse coupling strength (g syn ).Our results show that the spiral waves emerge for different combinations of R and g syn .We observe spiral waves when the value of the global phase order parameter is low (Z g < 0.7), the local phase order parameter value is high (Z L ≈ 1) and the amount of phase singularity is not greater than a maximum value (P max ≤ 20).We observe for low R values the desynchronized patterns.Moreover, there are many pattern transitions and bistable regions in which spiral and non-spiral waves are possible for higher R values.
The paper is organized as follows: In Section 2, we present our neuronal network model.In Section 3, we show the diagnostic tools.Section 4 presents our results.In the last Section, we draw our conclusions.

Neuronal Network Model
We consider a network composed of coupled adaptive exponential integrate-and-fire (aEIF) neurons [42].The dynamics of the neuronal network is given by where V i is the membrane potential, w i is the adaptation variable, g i is the synaptic conductance of the neuron i, I is the injected current, C is the membrane capacitance, g L is the leak conductance, E L is the resting potential, ∆ T is the slope factor, V T is the threshold potential, V rev is the potential reversal, τ w is the adaptation time constant, a is the level of subthreshold adaptation and τ g is the synaptic time constant.The parameter ∆ T controls the sharpness of the initial phase of the spike [45].This model improves the standard leaky integrate-and-fire by adding the adaptation current and the exponential term in the membrane potential.The adaptation current describes the slow activation and deactivation of some potassium ionic channels [46].Besides that, a spike threshold mechanism represented by the exponential term describes the fast arising when a action potential is generated [43,42,47].For ∆ T → 0, the neuron model becomes a standard leaky integrate-and-fire neuron model [48].If V i reaches the threshold V peak at a certain time t, the following reset conditions are applied: In our simulations, we consider the parameter values in Table 1.We constructed a two-dimensional excitatory neuronal network with 1000 µm X 1000 µm dimensions and N = 17, 324 neurons in which the distance between each neuron is 7 µm in x-axis and 8 µm in y-axis [41].The neuron i is connected with its neighbors that are within a radius R (M ij = 1), as displayed in Figure 1(a).

Diagnostic tools
To identify spiking and bursting firing patterns, we utilize the mean coefficient of variation (CV) of the neuronal inter-spike interval (ISI) [49,50], that is given by Table 1.Description and values of the parameters in the AEIF system used in the simulations [51].where σ ISI is the standard deviation of the ISI normalised by the mean ISI [52].The ISI is defined by ISI m = t m+1 − t m , where t m is the time of the m-th neuronal spike [53].Spike pattern produces CV < 0.5, meanwhile, burst pattern produces CV ≥ 0.5 [50].The mean firing rate is calculated by means of F = ISI −1 [50].

Parameter Description
To determine synchronous behaviour of the whole network, we consider the global phase order parameter, that is defined as [54] where φ j is the j-neuron phase, given by t m j corresponds to the time when a m-th spike (m = 0, 1, 2, . ..) of a neuron j happens (t m j < t < t m+1 j ).When the network is totally synchronized Z(t) = 1, total desynchronized states produce Z(t) = 0 and when 0 < Z(t) < 1 the network is partially synchronized [54,53].To better understanding of the network dynamics, we calculate the time-average of the order-parameter, that is given by where t fin − t ini is the time window.We use a time window equal to 5 s to compute Z to avoid transient effects.
We consider the local phase order parameter to study the local synchronization of each neuron.We separate the 1000 X 1000 µm network in 25 boxes in each dimension, totalizing 625 boxes, then each box contains a significant amount of neurons.The boxes are labeled by the subindex (i, j), where i indicates the x-axis and j the y-axis position of the box, i, j = [1,25].The size of each box is 40 × 40 µm and each one of them contains approximately 28 neurons.The local phase order parameter is given by [55] where k = √ −1, A is the total number of boxes, X i is the x-axis position of the box, and Y i is the y-axis position of the box.The mean value of Z i,j is the average of the local phase order parameter of the network is Figure 1(a) displays an illustrative scheme of the connections in our neuronal network, the neurons connect by means of chemical synapses.There are connections among neurons within the radius R, otherwise there is no connection.The radius origin is in the neuron position.Figure 1 (b-d) exhibit the time evolution of the membrane potential of three different neurons.In panel (b), for a low coupling strength g syn = 0.001 nS, the neurons are not synchronized.Increasing the coupling strength to g syn = 0.1 nS, the neurons fire in a synchronized spike pattern, as shown in panel (c).The panel (d) shows neurons synchronized and firing according to a burst pattern for g syn = 0.2 nS.The values of the local and global phase order parameters are exhibited in panel (e), where the circles, squares and triangles correspond to the firing patterns of the panels (b-d).The neurons in the network synchronize (locally and globally) for small values of g syn and the global order parameter varies much more than the local order parameter.The panel (f) displays the CV as a function of g syn .The CV value below the dashed line indicates the spike fire patterns.There is a maximum value of coupling strength before the transition from spike to burst behaviors, g syn = 0.2 nS.

Results
A bidirectional neuronal network is able to exhibit spiral wave activities [14] with different amounts of phase singularities (PS) and spiral waves.PS is the spiral core and is composed of asynchronous neurons.Each spiral wave rotates around its own PS and can share the same core.We observe that the local order parameter is below 0.7 in the PS.We define a Z i,j threshold to characterize the phase singularities.In this work, we focus on studying the correlation of spiral wave emergence with the network parameter R and g syn .We propose a reliable method to distinguish spiral wave from other travelling waves in the cortex.We use the diagnostic tools shown in the previous section.
Figure 2 displays three patterns of waves observable in the network and the local phase order parameter.The panels (a), (b) and (c) exhibit the phase (color bar) of each neuron in a time instant and the panels (d), (e) and (f) show the local order parameter of each wave pattern.The ring waves are shown in Figure 2(a).The wave pattern produces some dark orange dots in the panel (d), which are the origin of where the ring waves propagate from.The dark orange spots are not PS due to the fact that the local order parameter is not below 0.7.The spiral wave is displayed in the panel (b) and its centre is shown in the panel (e), by the black region.In the panel (e), we calculate the local order parameter of the spiral wave, where the black area corresponds to the phase singularity (PS).The mean local order parameter (Z L ) does not decrease significantly with the presence of the phase singularity, due to the high quantity of local synchronous neurons.For higher R values, it is possible to observe a burst spiral wave, shown in Figure 2(c).There are some border abnormalities in this case.In Figure 2(f), the border has low values of Z L , however, we do not count it as PS due the fact that there is no spiral wave rotating around it.To avoid mistakes in the PS counting, we discard the borders.High values of R increase the border effects due to the high connectivity of each neuron, although R influences the firing pattern.The panel (a) shows the parameter space of the global phase order parameter.The dark regions correspond to the asynchronous activity (Z g ∈ [0, 0.7]) and the brighter regions are the synchronized states of the whole network.The network shows synchronous behavior for high values of R. The value of g syn also plays an important role in the synchronization.Increasing it, we verify transitions from asynchronous to synchronous and vice-versa.The local phase order parameter is displayed in the panel (b).The coupling strength is of great importance for the emergence of synchronous states.The CV suffers changes with the increase of R and g syn .It is identified transitions from spike-to-burst, shown in the panel (c).The panel (d) exhibits the average firing rate, there is a quick increase for high values of R and g syn .In some regions of high g syn , it occurs the synchronization of bursting, which is strongly related to pathological cases, such as epilepsy seizures [50].The panels (a) and (b) are quite different in some regions, showing regions in which the network is globally desynchronized meanwhile local neurons are synchronized, Z g < 1 and Z L ≈ 1.This result characterizes waves in the network, however, it is not enough to distinguish the wave types.It is necessary more conditions to characterize spiral waves.The spiral waves have a particular characteristic, that is the phase singularities (the centre which the waves rotate around).However, high amounts of PS are related to asynchronous activity.We establish a maximum value of PS to avoid asynchronous activity.We propose the following conditions to detect spiral waves: (i) the mean global order parameter has to point out asynchronous activity, Z g ≤ 0.7, (ii) the mean local order parameter has to indicate a local-synchronized activity, Z L ≈ 1, (iii) there must be at least one phase singularity and PS is characterized by a low local phase order parameter, Z i,j ≤ 0.7, (iv) the number of PS is limited to a maximum number (P max ), 0 < PS ≤ P max .We consider P max = 20.
Conditions (i) and (ii) detect wave activity in the network, however, they do not describe the type of wave.The main difference among the wave types is their local phase order parameter pattern.For example, in the ring wave pattern, the local order parameter does not exhibit PS, as shown in Figure 2 (d).The presence of PS is characteristic of spiral waves, then at least one focus is necessary, condition (iii).Multi-phase singularities can be observed and the network becomes asynchronous when there are too many ones.We consider a maximum number of PS, condition (iv).These conditions all together indicate spiral wave in the network.The color scheme in Figs.The pattern parameter space of R × g syn is shown in Figure 5(a), where the colors yellow, gray, blue and red correspond to the spiral waves, non-spiral waves, asynchronous and synchronous patterns, respectively.For R ≤ 20 µm, it only generates asynchronous activity.Hence, the radius of connection has great influence in the type of pattern.There are some non-spiral waves regions in the range 40 ≤ R ≤ 35 µm and 0.25 ≤ g syn ≤ 0.05 nS, indicating a coexistence of synchronous and non-spiral waves pattern.It is important to highlight that the radius of connection is about 70 µm or greater in the human neocortex [41].The emergence of burst spiral is intrinsically correlated with long range connections (R ≥ 60 µm) and intense coupling strength (g syn ≥ 0.12 nS).The greenish and pinkish hexagons are the starting and ending point of the hysteresis analysis, respectively.Hysteresis diagrams of Z g and CV are displayed in Figs.5(b) and 5(c), respectively.Asynchronous and synchronous initial conditions are considered in the purple and green squares.The CV bistable region coincides with one of the regions in the Z g , both are observed in the g syn range [0.11, 0.13] nS.Two different outcomes are possible, that are burst synchronization for asynchronous initial conditions and spikes desynchronization for synchronous initial conditions.There is one more bistable region in the panel (b), in which the spike desynchronization is easier achieved for synchronous initial conditions.

Conclusions
In this work, we investigate the emergence of spiral wave patterns within a bidirectional network composed of pyramidal neuron cells.The individual dynamics of each neuron is described by the adaptive exponential integrate-and-fire model.The roles of spiral waves in cognitive activities is not fully understood.In this way, comprehending how network parameters influence the emergence of spirals and developing the capability to detect them is crucial for advancing our understanding of their role in cognitive activities.Hence, we investigate how network parameters influence the emergence of spirals.Additionally, we introduce a novel method capable of identifying chimera spiral waves.
We show that the network composed of coupled pyramidal neurons generates different spatial-temporal patterns.The local and global phase order parameters together have shown to be great tools of wave detection.We observe that wave activity generate a specific range of values for both mean order parameters.Locally the neurons are synchronized for almost all values of R and g syn , indicating that waves do not disrupt the local dynamics.The waves greatly affect the global dynamics and makes the global order parameter goes to values close to 0, characterizing desynchronization in the network.The low impact that waves have in local synchronization is due to the firing of the closest neighbors when the travelling wave arrives at the neuron.In the entire network, the passage of wave helps the local neurons to fire meanwhile the neurons that are not at the wave location do not receive a firing support.These characteristics hold great power to identify waves in the network, thus we suggest a diagnostic capable of identifying waves.They are identified by a low global and high local order parameters.Applying this criteria, we identify some waves.Mostly spike waves are found for radius lower than 55 µm.The type of wave is not pointed out by this diagnostic, then more criteria to identify spiral waves are necessary.
The identification of spiral waves can not be done by using only the condition proposed earlier, due to the fact that it does not specify the wave type.The spiral waves generate a specific abnormality in the phase producing very low local order parameter in it.The neurons localized in this abnormality are desynchronized and are the core of the spiral, called phase singularity (PS).The neurons in the PS exhibit a low value in the local order parameter inferior to 0.7, then we propose that PS are characterized by Z i,j < 0.7.It is not necessary only one PS to exist spiral waves, multi-phase singularies are possible.However, a great amount of PSs promote desynchronized activity.We define a maximum number of PSs in the network, P max = 20.We propose four conditions which must be met to identify spiral waves, these conditions regard the range of both phase order parameters and the amount of PS.The synchronized states are defined by high global and local order parameters, non-spiral waves are characterized by (Z L ≈ 1, Z g < 0.7 and PS = 0, spiral waves exhibit the same conditions of non-spiral waves with the addition of PS > 0 and desynchronized states are characterized by Z L ≈ 0 and Z g ≈ 0.
The spiral waves are generated by different combinations of R and g syn .This wave pattern is not only restricted for spike firing, but also for burst spiral.With regard to the burst spiral pattern, it is not as easily observed as spike spiral.The burst spiral appearance is restricted for specific values of R and coupling strength, for example R = 75 µm and g syn = 0.14 nS.The emergence of non-spiral wave is related to the initial conditions.To understand the initial conditions influence in the system, we analyze the hysteresis diagram for a constant radius and vary the g syn .Desynchronous behavior is easier to achieve if the initial conditions are desynchronized.We find a region of bistability in the range g syn =[0.11, 0.13] nS.Depending on the initial conditions, it is possible to verify the existence of spike or burst waves.The identification of bistable state is crucial due to its strong correlation with epileptic seizures [56].

Figure 1 .
Figure 1.The panel (a) displays a schematic representation of the two-dimensional network composed of coupled excitatory neurons.The panels (b-d) exhibit the time evolution of the membrane potential V for three different values of g syn .We observe (b) desynchronyzed spikes for g syn = 0.001 nS (blue circle), (c) synchronized spike for g syn = 0.1 nS (blue square), and (d) synchronized burst for g syn = 0.2 nS (blue triangle).The panel (e) shows the global and local Kuramoto order parameter as a function of g syn and the panel (f) exhibits CV as a function of g syn .In panels (b-f), we consider R = 70 µm.

Figure 2 .
Figure 2. (a) Spatial pattern of a synchronyzed spikes (R = 55 µm and Z g = 0.94), (b) spiral wave pattern (R = 64.5 µm and Z g = 0.15), (c) burst spiral wave pattern (R = 75 µm and Z g = 0.36), (d) local order parameter for (e) local order parameter of spike spiral wave and (f) local order parameter of burst spiral wave of each neuron i.We consider g syn = 0.14 nS.

Figure 3
Figure 3 exhibits the parameter space of the four diagnostic tools for R = 75 µm.The panel (a) shows the parameter space of the global phase order parameter.The dark regions correspond to the asynchronous activity (Z g ∈ [0, 0.7]) and the brighter regions are the synchronized states of the whole network.The network shows synchronous behavior for high values of R. The value of g syn also plays an important role in the synchronization.Increasing it, we verify transitions from asynchronous to synchronous and vice-versa.The local phase order parameter is displayed in the panel (b).The coupling strength is of great importance for the emergence of synchronous states.The CV suffers changes with the increase of R and g syn .It is identified transitions from spike-to-burst, shown in the panel (c).The panel (d) exhibits the average firing rate, there is a quick increase for high values of R and g syn .In some regions of high g syn , it occurs the synchronization of bursting, which is strongly related to pathological cases, such as epilepsy seizures[50].The panels (a) and (b) are quite different in some regions, showing regions in which the network is globally desynchronized meanwhile local neurons are synchronized, Z g < 1 and Z L ≈ 1.This result characterizes waves in the network,

Figure 3 .
Figure 3. (a-d) Parameter space of the network.The color code represents the values of (a) the global order parameter, (b) the local order parameter, (c) the mean coefficient of variation (CV) of the neuronal inter-spike interval and (d) the mean firing rate.

Figure 4 .
Figure 4.In (a) and (b) we plot, respectively, the global phase order parameter and the number of phase singularity as function of g syn for R = 75µm.The colors in (a) and (b) represent the dynamical pattern of the network.

Figure 5 .
Figure 5. (a) Parameter space of the network activity pattern, where each color represents one type of pattern.The yellow region corresponds to the spiral wave, the gray area is related to the non-spiral wave (for example ring pattern), the blue region denotes the desynchronized pattern and the red area corresponds to the synchronous pattern.The panels (b) and (c) show the hysteresis diagram for R = 75 µm, (b) global order parameter and (c) coefficient of variation.Asynchronous initial conditions are considered in the purple squares, meanwhile, synchronous initial conditions are used in the green squares.The greenish and pinkish hexagons exhibit the points which start and end the hysteresis cycle, respectively.