Investigation of Details in the Transition to Synchronization in Complex Networks by Using Recurrence Analysis

The study of synchronization in complex networks is useful for understanding a variety of systems, including neural systems. However, the properties of the transition to synchronization are still not well known. In this work, we analyze the details of the transition to synchronization in complex networks composed of bursting oscillators under small-world and scale-free topologies using recurrence quantification analysis, specifically the determinism. We demonstrate the existence of non-stationarity states in the transition region. In the small-world network, the transition region denounces the existence of two-state intermittency.


Introduction
Many natural phenomena can be modeled and studied through a mathematical approach. Especially, complex networks are useful for analyzing problems involving physical, biological, chemical, engineering, and even social perspectives [1,2]. In this way, coupled oscillators are able to investigate a large class of dynamical systems in a theoretical, computational or even experimental field [3][4][5][6]. As an example, neural networks can be understood as coupled oscillators where it is possible to associate the bursting neuron to a phase oscillator [7].
In the scenario of complex networks, the neural system can be modeled on two scales. In the internal connection scheme, each network node can be understood as a neuron and their connections as the edges [8], which are able to simulate a single network, as used in many works [9][10][11][12][13]. On the other hand, considering the inter-networks connection scheme, it is possible to consider a neural system composed of different sub-areas, so each sub-network can be understood as a node and their connections as the edges, building a network of networks [14][15][16][17][18].
The connection architecture of complex networks is very important in the dynamical properties observed at a global level of behavior. Regarding neural networks, many topologies are considered, such as small-world, scale-free, and random ones [9,16,19] where a transition from unsynchronized to synchronized states is observed. The role of connection architecture is very important to the paths to synchronization, where different phenomena can be observed as a function of topology [20][21][22]. In real neural systems, the characteristics of these topology schemes are observed [23][24][25], which motivates the investigation of the influence of topology on the dynamical properties of the networks.
Regarding complex networks, it is known that this kind of system can show emergent behavior, where the global behavior observed is richer than the sum of the individual element behaviors. In this way, the existence of non-monotonic transitions to synchronization as a function of coupling strength in neural networks [26][27][28][29][30], where non-stationary states can be noticed, has been reported in the literature. In some cases, in the transition region, on-off intermittency in the two states has been observed, where the network displays the existence of two locally stable states but globally unstable ones [18,26], as defined in [31]. In [9], the dynamical properties regarding synchronization and transition characteristics are studied as a function of the connection architecture, with both small-world and random topologies being considered. In the present paper, we extend the analysis and consider the scale-free topology, since there are topological differences between the connection schema. In the small-world network, there is a local connection structure plus a non-local one. On the other hand, in the random network, there is homogeneity in the connectivity distribution, which consists of an assortative network and composes a very different scenario in comparison to the scale-free network, which generally consists of hubs of connections and forms a disassortative network [32,33]. In fact, these differences may influence the dynamical properties of systems regarding synchronization, as observed in [20][21][22]33].
In order to simulate the neural behavior, we consider the coupled map developed by Rulkov [34]. Using this model, it is possible to reproduce bursting behavior, which is characterized by a sequence of chaotic spikes followed by a period of resting [35]. This kind of neural activity is observed in real neural systems, as reported in [36][37][38]. The building of networks involves two different topologies: small-world, obtained through the Newman-Watts route [39], and scale-free (power law distribution of connectivity), obtained through the Barabasi-Albert approach [40], since these topologies characteristics are observed in real neural systems.
To perform the numerical analyses of dynamical properties regarding phase synchronization and (non-)stationarity of the transition region, the Kuramoto order parameter [3] and recurrence quantification analysis (RQA) [41,42] are used. In RQA, the determinism is used, which evaluates the ratio of recurrent points that belong to diagonal structures. To use the Kuramoto order parameter, data from each neuron that composes the network is necessary since a phase is associated with the bursting behavior of all neurons. The use of recurrence quantification only requires a time series that characterizes the dynamical systems. In the case of networks, the determinism is evaluated from the mean field time series. It is known that the mean field of a phase synchronized network has a "periodic" oscillation where the amplitude is bigger than in an unsynhcronized case, as observed in [11,26]. This approach makes an experimental validation possible since the recurrence quantification analysis is able to analyze experimental time series [43][44][45][46].
In this paper, we investigate a network composed of 1024 chaotic bursting neurons simulated through the Rulkov map in terms of phase synchronization and transition region characteristics as a function of coupling strength. Here, we focus on the small-world and scale-free topologies and their influence on the dynamical properties of the neural networks. We show that the transition to phase synchronization depicts non-stationary characteristics; however, the transition occurs for smaller values of coupling strength in the scale-free network. On the other hand, the existence of two-state on-off intermittency is observed for the small-world network.
The paper is divided as follows. In Section 2, the Rulkov map and the bursting behavior are described. In Section 3, the connection architectures are shown, and the small-world and scale-free building approaches are described. In Section 4, the synchronization and intermittency quantifiers are presented. In Section 5, the results and a discussion are presented, which support the conclusions in Section 6.

Rulkov Map
The Rulkov map is able to reproduce different dynamical behaviors, such as chaotic spikes and a set of chaotic bursts [47]. The bi-dimensional Rulkov map can be described as where x i and y i are the fast and slow variables of the ith neuron. The set of parameters (α = 4.1, β = γ = 0.001) is chosen in order to obtain the bursting behavior, following [34]. The coupling term I i represents the influence of the other neurons in the ith neuron and is given by where ε is the coupling strength, χ is the normalization factor given by the average number of connections in the network [48], and a i,j is the adjacency matrix element where a ij = 1, if i and j are connected, 0, if i and j are not connected.
Here, we consider a small-world and a scale-free networks. Figure 1 depicts an example of the dynamical behavior of the fast and slow variables of the Rulkov map. Panel (a) represents x where bursting behavior can be observed, which is characterized by a sequence of chaotic spikes followed by a resting time. Panel (b) depicts y as a function of t, where each time that a burst starts, the slow variable assumes a maximum. For the interval of coupling strength (ε) used, the bursting behavior is maintained. Here, panel (a) depicts the fast variable of the system (x), which can be understood as the membrane potential of the neuron. Panel (b) depicts the slow variable of the system (y) and is useful for evaluating the oscillator phase since y assumes a maximum at every burst start.
Since the fast variable x reproduces the bursting behavior, it is possible to understand x as the neuron membrane potential, and the mean field potential can be described as [9] where N = 1024 is the number of oscillators (neurons) in the network.

Connection Architecture
The topology is very important to the dynamical properties of the networks. In fact, in neural networks, the connection scheme can be even more important than the structural variants of the neurons [49]. In this way, many results show that neural systems can display different topological properties, with small-world and scale-free characteristics being observed in real systems [50][51][52][53].
To build the scale-free network, we used the Barabasi-Albert approach [54], which can be described by two processes: • the network is expanded by the addition of new nodes; • the nodes added preferentially build connections to nodes that are already well-connected.
This approach is able to build a network where the distribution of the connectivity probability follows a power law: P(k) ∼ k −ν with ν > 1.
Regarding real neural systems, experimental evidence indicates that some brain activities may depict scale-free properties where the average number of connections is ≈ 4 [52,53]. Moreover, it is known that the human functional network may present scale-free characteristics where there is a power-law connectivity distribution following an exponent 2.0 ≤ ν ≤ 2.2 [52].
The scale-free network used in this paper has N = 1024 nodes with n = 4088 connections and ν = 2.20, which gives an average number of connections of χ sf = 4, which is used as the normalization factor in Equation (3).
In order to obtain the small-world connection matrix, the Newman-Watts route is used [39], where non-local connections are added in a second neighborhood regular network with a probability of p. The network is composed of N = 1024 nodes, and the total number of connections is given by where p = 1 leads to a globally connected network. Here, the number of local connection is 4096, and by using p = 5.3 × 10 −4 , 554 non-local connections are obtained, which leads to n = 4650 connections. For this network, the average number of connections is given by χ sw = 4.53.
It is possible to evaluate the average path length (L) and the clustering coefficient (C) [50] of the small-world network. For the matrix used in this paper, C sw = 0.3585 ∼ 10 −1 and L sw = 6.012 ∼ 1 are obtained using the NetworkX library [55]. However, it is possible to evaluate these quantifiers for an equivalent (with a similar number of connections) random network using the expressions L random ∼ ln (N)/ ln (n/N) and C random ∼ n/N 2 [56].
Defining the merit variable Γ = λ C /λ L where λ L = L sw /L random and λ C = C sw /C random [57,58] is possible to analyze the small-world existence condition. If Γ > 1, then the network considered has a small-world topology. For the network used here, C random ∼ 10 −3 and L random ∼ 1, which leads to Γ ∼ 10 2 , confirming that the networks used have a small-world topology. Figure 2 depicts the graph representation of the networks used in this paper. Panel (a) shows the small-world network, and panel (b) represents the scale-free one. One difference between the networks consists of the degree of connectivity, since in the small-world network the most connected neuron has 9 connections, while in the scale-free case this number is bigger than 140. Besides this, in the scale-free case is observed the formation of hubs characterized by larger connected neurons, which leads to non-homogeneity in the connectivity degree, while in the small-world network a different scenario is noticed, and the degree of neurons connectivity is similar.

Quantifiers
The synchronization of chaotic systems has been extensively studied since the last century [50,[59][60][61]. Here, we focus on the phase synchronization of networks composed of bursting neurons. There are different approaches to obtain the phase of chaotic oscillators; however, we have used the slow variable y to identify every time a burst starts. In this way, the phase is increased by 2π every maximum of y, and a continuous variation of the phase is obtained through [60,62] where t k,i is the time when the kth burst of the ith neuron starts. It is important to emphasize that the bursting behavior is maintained for the entire interval of coupling strength used (ε), which means that phase can be obtained from the time series of the maximum of y without reconstruction of the dynamics in the higher dimensional phase space [60].
To evaluate phase synchronization based on phase θ, the Kuramoto order parameter is used [3]. Considering the phase of all neurons, it is possible to define If the system is (not) in a phase synchronized state, then R → 1 (R → 0).
To obtain the synchronization level as a function of the coupling strength, the time average of R(t) can be computed for each value of ε. The mean value of the Kuramoto order parameter is given by where t f is the total simulation time, and t 0 is the transient time.
Recurrence quantification analyses are useful for analyzing the synchronization characteristics of dynamical systems [9,26,42,63,64]. The original work about recurrence analyses is based on a visual approach to identifying dynamical properties [65]. On the other hand, it is known that the mathematical structure of recurrence plots (RP) offers information about the dynamics of a system [41,42].
Here, we perform recurrence analyses based on the time series of the mean field of the network, described by Equation (4). In this case, the recurrence matrix is defined as where Θ is the Heaviside function, δ is the recurrence threshold, and S is the size of time series analyzed. Based on a time series of size S, of the mean field potential of the network described by Equation (4), it is possible to obtain the recurrence matrix through the use of Equation (9). The procedure consists in analyzing whether each point of the time series is close enough to another point (the difference between them must be smaller or equal to δ). After considering the entire time series, a matrix (S × S) is obtained where if a state x a is (not) recurrent to another x b , R ab assumes one (zero). The visual approach to analyzing the system consists of relating the recurrent (not recurrent) point to a black (white) dot [65]. In order to fix δ ∈ [0, 1], the time series of the mean field must be normalized.
In order to analyze the synchronization of the network, the diagonal lines in the recurrence matrix are very important since their distributions are related to the regularity of the trajectories. A diagonal line of length is understood as a segment of the trajectory rather close to another segment of the trajectory during time steps in a different time [41].
In this way, the better quantifier for the synchronization analysis is the determinism, which gives the ratio of recurrent points that belong to diagonal lines over all recurrent points in the recurrence matrix defined by Equation (9). This method is useful since the mean field of a phase synchronized network, depicted by Equation (4), has "periodic" oscillations, as observed in [26,30].The determinism is defined as where min is the minimum length to consider a diagonal line. P( , δ) is the probability distribution function (PDF) of the diagonal lines.
In a similar approach that uses the Kuramoto order parameter, it is possible to evaluate the mean value of the determinism using

Results and Discussions
The numerical results have been obtained from initial conditions that follow a random distribution in the intervals described by the fast and slow variables (see Figure 1), which avoid any initial trend. The transient time is given by t 0 = 100,000.
An important point to note is the dependence of the recurrence quantification analyses (RQA) on the recurrence parameters, an particularly, the recurrence threshold (δ) and minimum diagonal length ( min ). If δ → 1, all points will be considered recurrent, and the determinism becomes saturated [41]. Similar behavior is obtained if min → 1, since all points will be considered as diagonal structures [41]. In order to use RQA to investigate dynamical properties regarding phase synchronization and/or non-stationarity, it is necessary to optimize the quantifiers through the choice of recurrence parameters [9,43]. The choice of recurrence threshold parameter follows [43],where the condition is described by d[∆(δ)]/dδ = "a maximum" leading to δ = 0.11, which results in the highest sensitivity of the quantifier determinism(∆). The minimum diagonal length ( min = 35) is chosen in order avoid small diagonals, and the determinism (∆) is better able to distinguish the phase synchronized states [9]. The determinism is evaluated using a moving window of 10,000 points.
The phase synchronization main scenario is depicted in Figure 3, where the synchronization characteristics are studied as a function of the coupling strength. Here, the total time of simulation is given by t f = 250,000. Panels (a) and (b) depict the mean value of the Kuramoto order parameter ( R ) and the mean value of the determinism ( ∆ ) for the scale-free network. Panels (c) and (d) depict the same analysis for the small-world network. Here, 20 different seeds for the initialization of the system are considered, and the vertical magenta bars indicate the dispersion (standard deviation) over initial conditions. For both networks, a clear transition from unsynchronized to phase synchronized state is observed; however, in the scale-free network, the transition occurs for smaller values of coupling strength in comparison to the small-world network. The dispersion over initial conditions is bigger at the transition region for both cases; however, for the small-world network, the phenomenon is more visible. A similar scenario is observed in [9]. In the scale-free network, the transition occurs where 0.0025 < ε < 0.012, where the increase in the coupling strength makes the system reach the phase synchronized state. In the small-world network, the transition region is characterized by 0.005 < ε < 0.017, and phase synchronization is observed for higher values of ε. This fact indicates that the scale-free network reaches the phase synchronized states for smaller values of ε in comparison to the small-world network. A similar scenario is observed in [66], where small-world networks reach synchronous behavior slower than other topologies. However, it is important to emphasize that the coupling term given by Equation (3) depends on the coupling strength (ε) and normalization factor χ. In this case, the normalization factor for each network is given by χ sf = 4 (scale-free) and χ sw = 4.53 (small-world), which leads to a relative factor of 1.13 between them. On the other hand, the difference in the coupling strength regime for the occurrence of the transition and phase synchronization is bigger than 1.13, so it is possible to conclude that the scale-free network reaches the phase synchronized state for a smaller coupling strength than the small-world one. Considering the small-world and random networks, it is possible to observe a higher level of synchronization in the random case [9,16]. Here, a similar scenario to the small-world and scale-free case is obtained; however, a new phenomenon regarding the critical value of ε is observed, as previously mentioned. In fact, the scale-free network presents a smaller value of ε to reach synchronization in comparison to the small-world and random cases, as studied in [9]. A similar scenario is observed in [20].
In order to investigate the details regarding synchronization characteristics, Figure 4 depicts the fast variable x for all neurons in the network as a function of t. Panel (a) represents the scale-free network with ε = 0.002, panel (b) ε = 0.005, and panel (c) ε = 0.040. Panel (d) represents the small-world network with ε = 0.002, panel (e) ε = 0.011, and panel (f) ε = 0.040. In the unsynchronized states (panels (a) and (d)), the neurons start their bursts without any coherence; however, in the transition region (panels (b) and (e)), the formation of horizontal structures is noticed. Finally, in panels (c) and (f), the phase synchronization where the horizontal structures denounce the spatial-temporal coherence of their bursts can be noticed. The mean field of the networks (Equation (4)) is depicted in Figure 5. It is important to note that the mean field consists of more easily experimentally accessible data than the individual signal. As observed in [12,30,67], the amplitude of the mean field of the network is related to the phase synchronization. Here, the scale-free (panels (a), (b), and (c)) and small-world networks (panels (d), (e), and (f)) are considered. In the unsynchronized cases (panels (a) and (d)) where ε = 0.002, the amplitude of the mean field V is vanishing. However, increases in the coupling strength (panel (b)-ε = 0.005 and panel (e)-ε = 0.011) lead the network to the transition region, where a small amplitude of V is observed. Importantly, the amplitude of the mean field varies as a function of t, as observed in [26], which indicates the possible existence of intermittency. Finally, for higher values of coupling ε = 0.040 (panels (c), and (f)), the mean field depicts an oscillatory behavior with a higher amplitude, and the frequencies are related to the bursting activity, which indicates phase synchronization behavior.
The "periodic" behavior of the mean field time series makes possible the use of RQA to evaluate the synchronization characteristics of the networks [9,26,67]. As depicted in Figure 3, the mean value of the determinism ( ∆ ) is able to distinguish the level of phase synchronization. Figure 6 depicts the recurrence plot obtained from the recurrence matrix described by Equation (9), where the black dots are related to the recurrent points and the white dots to the non-recurrent ones. Here, the recurrence matrix is obtained from the mean field time series described by Equation (4) and depicted in Figure 5. Panel (a) of Figure 6 depicts the recurrence plot of an unsynchronized state, while panel (b) depicts the recurrence plot of a phase synchronized one. It is noticed that in the phase synchronized case, the diagonal structures are much clearer, which leads to a higher value of determinism (∆) and makes its use possible for analyzing the synchronization characteristics of networks.  The determinism is able to distinguish between more (less) synchronized states with just the use of the mean field time series [26]. In this way, it is possible to analyze the time series of the determinism in order to obtain information about synchronization characteristics as a function of t. Following [9], the temporal standard deviation of the determinism time series for different values of coupling strength (ε) offers information about the intermittency between states with different synchronization levels.
where A represents the Kuramoto order parameter (R) or the determinism (∆), and T is the length of the time series considered. Here, a stationary time series of R or ∆ leads to a vanishing value of σ.
On the other hand, if the time series of R or ∆ display an intermittent behavior, where the network assumes different synchronization level states as a function of t, then σ assumes higher values [18]. Figure 7 depicts the temporal standard deviation of the Kuramoto order parameter and the determinism time series as a function of the coupling strength (ε). Here, panels (a) and (b) depict σ(R) and σ(∆) for the scale-free network, while panels (c) and (d) depict σ(R) and σ(∆) for the small-world network, where σ is normalized by its maximum value. For all cases, it is observed that increases in the coupling strength lead to σ increases, until a maximum. After reaching the maximum, increases in ε lead to σ decreases until it reaches a vanishing value. Similar behavior is observed in [9] , where the maximum of σ is located at the transition region of the coupling value.Here, the observed behavior is the same since the transition region depicted in Figure 3 is related to the region where σ depicts higher values. This is indicative of intermittent behavior [18]. Again, for the scale-free network, the transition region, and consequently the intermittent region, occur at smaller values of coupling strength in comparison to the small-world network (see discussion about Figure 3). In particular, the maximum occurs at ε ≈ 0.0045 in the scale-free network, and for the small-world network, the maximum is observed at ε ≈ 0.010. Figure 7. Normalized temporal standard deviation (σ), described by Equation (12), as a function of coupling strength for the scale-free network (panels (a) and (b)) and the small-world network (panels (c) and (d)). Here the time series of the Kuramoto order parameter (panels (a) and (c)) and the determinism time series (panels (b) and (d)) are considered. A higher value of σ indicates an intermittent behavior. The local maximum is observed in the transition region, indicating non-stationary transitions.
The determinism time series (∆) is depicted in Figure 8 in order to obtain more details regarding the network dynamical properties in the different situations, as depicted in Figures 3-5. As previously mentioned, the recurrence matrix (Equation (9)) is evaluated from the mean field time series. In this way, the diagonal structures observed in Figure 6b are related to trajectories that are close enough (smaller or equal to δ) at different times. Through the optimization of δ following [43] , it is possible to use the determinism to quantify the synchronization level of the network as a function of time. In Figure 8, scale-free (panels (a), (b), (c), and (d)) and small-world networks (panels (e), (f), (g), and (h)) are considered, where panels (a) and (e) represent the unsynchronized states with ε = 0.002. It is observed that the determinism assumes small values, which correspond to a unsynchronized states [9]. Panels (b) and (c) represent the states at the transition region for the scale-free network (ε = 0.0045 and ε = 0.0050), where it is observed that the determinism shifts between different values. This fact indicates that the network assumes different synchronization states as a function of t, characterizing an intermittent behavior. Similarly, panels (f) and (g) represent the states at the transition region of the small-world network (ε = 0.0105 and ε = 0.0110). Again, an intermittent behavior is observed, where the determinism shifts between different values; however, in this case, two plateaus of synchronization are observed, indicating the possible existence of two-state on-off synchronization [18,26,68]. In [9], a similar system is considered under small-world and random topologies. Despite the random case not showing the intermittence between two states, it is possible to observe differences in comparison to the scale-free case, since in the scale-free case, the existence of the higher state is less pronounced. Finally, panels (d) and (g) represent phase synchronized states for the scale-free and small-world networks (ε = 0.040), respectively,where a stationary signal is observed for high values of ∆ in both cases.  Figure 9 depicts the probability distribution functions (PDFs) of the determinism time series with a length of 10 7 points, since long-term data is necessary to evaluate the temporal behavior of the networks [26]. The first row, panels (a), (b), (c), and (d) represent the scale-free network while the second row, panels (e), (f), (g), and (h) represent the small-world one. Panels (a) and (e) represent the unsynchronized state of scale-free and small-world networks, respectively, where ε = 0.002. In this case, a uni-modal distribution with a small dispersion over the mean value of ∆ is observed, indicating the unsynchronized state [26].  Figure 8. The panels follow the same scheme as in Figure 8, where the first row refers to the scale-free network, and the second row refers to the small-world network. In the transition region, the distribution is not symmetric and indicates a non-stationary behavior; however, only in the small-world case is the existence of two states observed.
If the coupling strength is increased, the networks assume the transition states, where the temporal standard deviations (σ) of the determinism time series are higher for both topologies considered, thus an intermittent behavior is expected. In this way, panels (b), (c), (f), and (g) represent the states at the transition region for the scale-free and small-world networks. Here, in panel (b) ε = 0.0045 and in panel (c) ε = 0.0050, a non-symmetric distribution of ∆(t) is observed, which indicates a non-stationary situation. In panels (f) and (g), the PDFs of ∆(t) at the transition state of the small-world network is depicted, where ε = 0.0105 and ε = 0.0110, respectively. Here, a non-stationary behavior is again observed since there is a non-symmetric distribution. However, for these cases, the existence of two states with similar occurrence probabilities is observed in the presence of two peaks in the PDFs. A similar scenario is observed in [9,18,26]. In fact, in [9] , the distribution of ∆ with two peaks is observed in the small-world case. In the random case, a main peak and a shoulder in the distribution of ∆ is observed. A comparison between the random and scale-free topologies makes it possible to notice that in both cases, there is no two-state intermittence; however, the random case presents a more unsymmetric PDF of ∆, which indicates a higher level of non-stationarity.
Finally, panels (d) and (h) of Figure 9 are representative of phase synchronized global stable states for scale-free and small-world networks where ε = 0.040. The PDFs are uni-modal with a very small dispersion over the mean value, which indicates that the state is stationary and globally stable [26]. It is important to emphasize that for the scale-free case (panel (d)), the mean value of the PDF of ∆ is higher than in the small-world case (panel (g)), which suggest a higher level of phase synchronization.

Conclusions
The characteristics of synchronization in complex networks are very interesting and can be used for application in many nature phenomena [1,2,42,60]. Recently, studies have identified a non-stationary transition in bursting neurons networks [11,18,26] where, in some cases, the existence of two-state on-off intermittency is observed.
In order to analyze neural networks, we have used the Kuramoto order parameter [3], which uses the individual signal of all neurons in the network through the association of a phase based on the bursting activity. Moreover, we used recurrence quantification analysis [41,42,65], which is based on a time series that characterize the dynamical system. To analyze the mean field time series, we used the determinism, which is the ratio of recurrent points that belong to a diagonal structure.
Here, we have focused on a scale-free network composed of 1024 bursting neurons (oscillators) with 4088 connections and a small-world network composed of 1024 bursting neurons coupled trough 4096 local connections and 554 non-local ones. We noted a clear transition from unsynchronized to phase synchronized states where non-stationary behavior was observed in both cases. A similar scenario has been observed in [9,11,26,27,69]. Despite the higher number of connections in the small-world case, the transition region and the phase synchronized states were observed for smaller values of coupling strength in the scale-free case.
We have based our analyses in methodology proposed in [9], but here, we extended the study to a scale-free connection architecture since topology plays an important role in the dynamics of systems. It is important to emphasize that there are significant differences between small-world, random, and scale-free topologies. The small-world connection architecture is characterized by a regular local connections scheme with the addition of non-local ones. The random connection architecture is characterized by a more homogeneity in the connectivity, while the scale-free one is the opposite, since there are hubs of connections. The results have shown that the small-world case demonstrates the existence of two-state intermittence, as explored in [68]. Regarding the scale-free network, despite the non-stationary transition, this phenomenon is not observed. Similar results are obtained in [9] for random networks; however, it is also possible to observe differences between scale-free and random networks, as the random ones present a more non-stationary transition.
The study of complex networks has many application fields and is of great theoretical interest. Here, we have demonstrated that small-world and scale-free networks composed of bursting neurons can show non-stationary behavior at the transition to phase synchronization. These kinds of topologies find support in real neural systems, besides a large class of dynamical systems [1,23,50,52]. Thus, the present paper offers a theoretical/computational approach to understanding neural phenomena, since synchronization and intermittent behavior may be related to neural diseases [70,71] such as Parkinson's disease, autism, and Alzheimer's [72][73][74][75].