Synchronization malleability in neural networks under a distance-dependent coupling

We investigate the synchronization features of a network of spiking neurons under a distance-dependent coupling following a power-law model. The interplay between topology and coupling strength leads to the existence of different spatiotemporal patterns, corresponding to either non-synchronized or phase-synchronized states. Particularly interesting is what we call synchronization malleability, in which the system depicts significantly different phase synchronization degrees for the same parameters as a consequence of a different ordering of neural inputs. We analyze the functional connectivity of the network by calculating the mutual information between neuronal spike trains, allowing us to characterize the structures of synchronization in the network. We show that these structures are dependent on the ordering of the inputs for the parameter regions where the network presents synchronization malleability and we suggest that this is due to a balance between local and global effects.


I. INTRODUCTION
Modeling of natural phenomena through the use of coupled networks finds applications in various scientific areas [1,2]. Examples can be found in statistical physics [3], power grid distributions [4], ecology [5], and neuroscience [6]. This approach has revealed a variety of synchronization phenomena in coupled networks, such as complete synchronization, and phase synchronization (PS) [7][8][9]. Building on these findings, recent investigations have shown cases where the network susceptibility plays an important role so that small changes in the system can produce significant consequences [10]. Examples are in synchronization vulnerability, where small perturbations in nodes can lead to a desynchronization process [11], and also changes in a small number of initial conditions can lead to different dynamical features related to chimera states [12].
In neuroscience, synchronization has been found to be of fundamental importance: it has been observed in healthy behaviors, like memory [13], conscious processes [14], visual-motor behavior [15] and perception phenomena [16]. Also, either excess or lack of synchronization have been related to unhealthy behaviors, like seizures, generating epileptic episodes [17], Parkinson's disease [18], and autism [19].
One way to study synchronization is by the modeling of a neuronal network. There are several models to reproduce the dynamics observed in neurons [20]. A suitable model is the two-dimensional map developed by Chialvo [21], which mimics spiking behavior. Networks composed of neurons simulated by this model can present a diversity of dynamical phenomena, such as bistability with explosive synchronization (in small-world networks) [22] or even the emergence of synchronization patterns like clustering synchronization and anti-phase clusters, occurring due to the interplay of interacting sub-populations [23]. Besides the neuronal model, the connection topology is also an important factor for synchronization phenomena [24]. In general, it has been found that connection architectures composed of local structures do not facilitate higher degrees of synchronization [25]. However, long-range and global connection schemes do facilitate synchronization processes [26]. To analyze these phenomena, a coupling architecture given by a distancedependent power-law scheme can be used. In this kind of coupling, all neurons are connected, but not necessarily all connections are effective since the contribution of more distant nodes decreases with the increase of the power-law exponent. In fact, by varying this parameter, there is a continuous transition from global effectiveness (all neurons contributing equally) to local effectiveness (only first neighbors contributing). This topology has been studied in several contexts since the long-range interaction is observed in fundamental laws of physics, in coupled-oscillators networks [27], in biological networks [5], and in the connectivity of neurons [28].
In this paper, we focus on the synchronization properties of networks of Chialvo neurons coupled through the distance-dependent power-law scheme. We use the Kuramoto order parameter [29] to measure PS, with the neuronal spiking activities being associated with geometric phases. The main results consist of the description and analysis of the phenomenon we call synchronization malleability, in which the PS behavior of the network drastically changes as a consequence of variations of the neuronal inputs. We consider a sequence of neuronal input values with a uniform distribution and show that a shuffling process over this sequence can change the network from a highly phase-synchronized state to a highly desynchronized one. Besides, we show that these systems present diverse synchronization patterns due to the interplay between the coupling strength and the powerlaw exponent. For networks where the dependence on distance is weaker (closer to the global case), we observe a traditional transition from non-synchronized to phase-synchronized states as the coupling strength is increased. On the other hand, as we make the distantdependence stronger (closer to the local case), there ap-pear states where only parts of the network are phasesynchronized and there appear also diagonal spatiotemporal structures, or zig-zag states [30,31]. For sufficiently strong distance-dependence, these states start to dominate and the network no longer transitions to phase synchronization.
We also measure the mutual information [32] shared between neurons, evaluated through their spike-count [33,34], which is a measure of their PS degree [35]. With this, we obtain the functional connectivity of the network, allowing us to characterize its synchronization structures. When the information is effectively shared between neurons at both local (neighborhood) and global (long-range) levels, network PS is reached. Otherwise, in the cases where the network is non-phase-synchronized the sharing of information is effective only at the local level. We see that different shuffles of the neuronal inputs may facilitate the sharing of information at the global level, leading to network PS. However, they may also hinder it, resulting in a PS degree smaller than observed in the uncoupled case.
The paper is organized as follows: in section II, the local dynamics and the connection architecture are presented; in section III, the methodology is exposed; in section IV, the main results are shown and the discussions are presented to support the conclusions in section V.

II. THE MODEL AND CONNECTION ARCHITECTURE
To simulate the local dynamics, we consider the neuronal model [21] x where x i,t and y i,t are the activation and recovery variables of the i-th neuron, with the variable x i,t mimicking the membrane potential. K i is the input signal in each neuron, which acts as an additive perturbation [21], affecting the neuron's firing rate. The constants a, b, and c are parameters that control the dynamical behavior of the model, set to obtain the spiking behavior (a = 0.89, b = 0.6, and c = 0.28) [21]. I i,t is the coupling term between neurons and follows a distance-dependent powerlaw scheme described as where ε is the coupling strength, N is the network size, N = (N − 1)/2, α is the power-law exponent, or locality parameter, and η is the normalization factor given by It is important to notice that α controls the range of effective connections, which contribute significantly to the coupling term I i,t . When α = 0, a typical global coupling scheme is obtained, where the i-th neuron is effectively coupled with all the other N − 1 neurons. When α is increased, the distance-dependence becomes stronger, and more distant neurons contribute less. In the extreme case of α → ∞, the coupling scheme is characterized by a first-neighborhood topology. Figure 1 depicts the dynamics of one isolated neuron. Panel (a) shows x i,t as a function of t, characterizing a depolarization and repolarization process (spiking dynamics). Panel (b) depicts the recovery variable dynamics y i,t . For the entire paper, similar dynamical features are observed for all coupled neurons. The black circles indicate when the activation variable reaches the condition x i = 0.5 with positive first derivative, leading to a spike event. The sequence of spike times for one neuron is called its spike train, which we use in the evaluation of phases and production of raster plots. The input values K i are initially built to produce a uniform distribution where σ = 0.0035 is the coefficient of neuronal dissimilitude. Different values of σ were tested and similar results were obtained, as long as σ is kept sufficiently small to guarantee the spiking behavior [21]. For the simulations, we perform a random shuffling process over the K i values obtained from Eq. (5). This leads to a different ordering in the input values but keeps the distribution still uniform. We consider 30 different shuffled sequences, labeled as shuffling #1, shuffling #2, ..., shuffling #30. The relevant codes and sequences of input values ({K i }) can be found in the repository [36]. Similar results are also observed when K i are obtained from random generators considering the limits [0.03, 0.03 + σ]. Figure 2 shows the values of the inputs K i as a function of the neurons' index i. The black dots represent the case constructed following Eq. (5). The pink up triangles and the gray stars correspond to shuffling #1 and #2, respectively.

A. Kuramoto order parameter
The Kuramoto order parameter [29] is used to quantify phase synchronization (PS) between oscillators. To do that, we first define a phase θ i for the i-th neuron such that θ i increases by 2π for every spike. A continuous variation of θ i can be obtained through a linear interpolation [37] θ i (t) = 2πn i + 2π t − t n,i t n+1,i − t n,i , (t n,i ≤ t < t n+1,i ), (6) where t n,i is the time when the n-th spike occurs in the i-th neuron. The degree of PS of the network for each time t is, then, with i in this equation denoting the imaginary unit. The quantifier R ranges between 0 and 1. If R = 1 there is complete phase synchronization. Otherwise, if R → 0 the system can be non-synchronized, meaning spatiotemporal incoherence, or it can even display antiphase-synchronization. We take the time average to obtain the mean Kuramoto order parameter R : in which t 0 is the transient time and t f is the whole simulation duration. The Kuramoto order parameter can be evaluated for only a group of neurons in the network [38]. In this case, the network is divided into communities [38] and the mean Kuramoto order parameter computed for each community is given by where the phases θ k (t) are obtained from Eq. (6), N local = N/M is the number of neurons in each community when the entire network is divided into M communities, and Ω j is the set of neurons belonging to the j-th community. The partitioning is done such that the j-th community Ω j contains neurons with indices in the interval [jN local , (j + 1)N local ), with j = 0, 1, · · · , M − 1. Eq. (9) allows us to compare the PS in the communities and the entire network by defining where R local = 1 /M M j=1 R j is the average over all M communities. If δR → 1, the communities are phase synchronized but the entire network is not. On the other hand, if δR → 0, the communities and the entire network have similar dynamics, either non-synchronized or phase synchronized.

B. Spiking frequency
Given the phases θ i (t) associated to each neuron in the network, we evaluate the angular frequency related to the spiking activity as: In this way, to quantify the degree of frequency synchronization in the network, we evaluate the standard deviation of the spiking frequency ω i over all neurons in the network and then divide by its mean value: where ω = 1 /N N i=1 ω i is the mean frequency. In this sense, if κ(ω) → 0, neurons have the same spiking frequency and the network is in a state with frequency synchronization.

C. Information analysis
In this subsection, we describe the quantifiers we used to measure the information contained in the neuronal spike trains. The coding mechanism we consider is the spike-count code, according to which information is coded by the neuron in the number of spikes over some relevant time window [39]. In this case, we partition the spike train into bins of size ∆t and count the number of spikes per bin [33,34]. The choice of the parameter ∆t will be discussed later. Then, the information in the spike train of neuron N is calculated as its Shannon entropy [40]: where n is the number of spikes in each bin and p(n) is the probability of observing n spikes, estimated as the frequency of bins with n spikes in the time series. For two neurons M and N , we can calculate the amount of information shared by these two (or, the amount of information that one contains about the other) using the mutual information [32,41]: where p(m) and p(n) are marginal probabilities and p(m, n) is the joint probability of the values m and n. For a wide range of ∆t, an oscillating pattern in the mutual information between neurons is observed. We choose the ∆t that corresponds to the first maximum, following the principle of maximization of entropy [33] and giving us the clearest results. As a consequence, the bin size ∆t isn't necessarily the same for every parameter value.
The challenge with using M I is the estimation of the probability distributions. We use the method known as the plug-in, or direct method [43,44], which results in positively biased estimates. However, this positive bias can be reduced by increasing the number of samples [43,45]. We consider 500000 times for the analyses, which give us stable results: increasing this time did not alter results considerably.
At last, we collect the mutual information between every pair of neurons (M, N ) in a matrix (MI) MN = M I(M, N ) and consider this as the functional connectivity (FC) of the network.

IV. RESULTS AND DISCUSSION
We study a network composed of N = 525 spiking neurons following Eqs. 0.01, the network is non-phasesynchronized (dark blue tones). However, for α 1.5, the increase of coupling strength ε leads the network to phase-synchronized states (red tones), a traditional scenario observed in several contexts [24,29,37]. Increasing the locality parameter to the range 1.5 α 2.0, the transition to PS is maintained, but on average with a reduced degree (smaller R ) and a higher ε value being required to reach phase-synchronized states. At last, for α 2.0, the network no longer transitions to PS, since, for the entire interval of ε, R depicts blue tones. In this case, the increase of the coupling strength (ε 0.01) actually reduces the average degree of PS and the network depicts values of R smaller than the uncoupled case. Hence, we observe two distinct transitions between non-phase-synchronized states and phase-synchronized states: one induced by the increase of coupling strength and another by the decrease of α. In the former, an increase of ε leads to PS. In the latter, an increase of α terminates PS, as the network becomes locally coupled. A similar scenario was found considering different local dynamics [46][47][48].
An interesting behavior is observed in panel (b) of Fig.  3, where the dispersion of the networks' PS degrees is shown. In this case, higher values of χ(R) indicate that the system can depict different degrees of PS in different simulations. This is most pronounced in 1.6 α 2.3 and ε 0.020, where there is the transition to PS induced by a decrease in α (as seen in panel (a)). This region is called malleable, since the shuffling process in the input values K i leads to networks presenting different degrees of PS, and, consequently, different dynamical states for the same set of parameter α and ε. For other regions in the parameter space (ε 0.020 or α 2.3 or α 1.6), the quantifier χ(R) depicts small values, indicating that the network's dynamics is similar over simulations. We note that the synchronization malleability phenomenon occurs for intermediate values of α, between the extreme cases of local and global effectiveness. A similar phenomenon was observed in a small-world network of bursting neurons when some local connections are changed by non-local ones [49].
Panel (c) depicts the quantifier δR , given by Eq. (9). For ε 0.020 and 1.5 α 2.7 higher values are observed, indicating that there is PS in the communities, but not in the entire network. Particularly interesting is the case for 1.6 α 2.0, since the increase of ε, for a fixed value of α, makes δR transition from small (blue tones) to high values (red tones) and then back to small values (blue tones). For these regions, the same increase of ε makes the network transition from non-phase-synchronized states to phasesynchronized ones (see panel (a)), where δR indicates that this transition occurs with the existence of PS at the community level before the entire network reaches phase-synchronized states. For α 1.5, δR always depicts small values, and the transition to PS induced by the coupling strength does not pass through states with community PS. For high values of α, α 2.9, the quantifier shows also small values, since the network does not reach PS. These analyses were performed considering M = 15 communities, but similar results are obtained considering M = 7 and M = 5 communities. The absolute value changes, but the region in the parameter space where the quantifier depicts higher values is the same, which corroborates the main idea of the existence of phase-synchronization on a local level, but not on a global one. Details of the dynamical states can be observed in the raster plots (RP) of the network (see Figs. 5 and 6). Fig. 3 shows the quantifier κ(ω), described by Eq. (12). For the entire interval of α and ε 0.010, smaller values indicate that the network is not frequency synchronized. The increase of the coupling strength beyond ε 0.010 then leads the network to frequency-synchronized states. This phenomenon is observed for small values of α, where the network is on a phase-synchronized state (high value of R ), which is expected. However, for higher values of α, the network is also frequency synchronized, but not phase synchronized.

At last, panel (d) of
In order to analyze the details of the network's dynamics, Fig. 4 (a) depicts the maximum and minimum values of R , represented by the extremes of the filled area, for the 30 simulations as a function of the coupling strength ε. The same representation is used in panel (b) for the maximum and minimum values of δR . For α = 1.0 (blue lines with circles), an increase in ε makes the network reach PS, as R increases and δR is small for all simulations. This is expected since in Fig. 3

(b) χ(R)
shows a small value in this region.
On the other hand, for α = 1.8 (orange lines with squares), the network may assume significantly different synchronization states, as seen in the difference between the extremes of R . For 0.020 ε 0.080, with fixed parameters α and ε, the network can either be phase-synchronized, with R ≈ 0.90, or even be nonphase-synchronized, with R ≈ 0.05. The difference is simply due to different shuffling over the sequences of input values K i . This region also depicts a great difference between the maximum and minimum values of δR . Higher values of this quantifier occur with smaller values of R , in which case the network is only phase synchronized at the level of communities, but not globally. Otherwise, smaller values of δR occur with higher R , in which case the entire network reaches PS. This phenomenon, called synchronization malleability, is observed in the region where χ(R) has high values ( Fig. 3  (b)).
For α = 2.5 (green lines with diamonds), the network is always non-phase-synchronized, with R < 0.4. Furthermore, we note that some of the values of R for the coupled network are smaller than in the uncoupled case (ε = 0.0). In fact, in these cases the values of R are similar to the expected value for randomly distributed phases (R random ∼ 1/ √ N = 0.0436, considering N = 525) [9]. Besides, the increase of ε makes the maxima and minima of δR also increase, indicating that the communities of the network may be phase synchronized, even though the network as a whole may not.
Panel (c) of Fig. 4 depicts the Kuramoto order parameter R as a function of ε for the non-shuffled case, which is the sequence of input values K i given by Eq. (5) and represented by the black dots in Fig. 2. The red line (big diamonds) represents the case of α = 1.0, where the increase of ε makes the network transition from nonsynchronized states to phase-synchronized ones. Comparing with the shuffled cases (panel (a)), the nonshuffled case requires higher coupling strength to phase synchronize and even then it does so at a lower degree of R . Increasing the locality parameter to α = 1.8 (purple line with hexagons), the network is always nonphase-synchronized, a drastically different scenario than the shuffled cases, where the network can depict PS. At last, for α = 2.5 (brown line with crosses), the network also depicts small values of R and does not reach PS. We therefore see that the shuffling may facilitate the PS of the network or even hinder it. This limit case, represented by the non-shuffled values of K i , is just an artificial construction used to illustrate the effects of the shuffling process.
In panel (d) of Fig. 4, the dependence of the network to initial conditions is studied. It depicts the R of the network for two representative shuffled input sequences {K i } shown in Fig. 2, labeled as shuffling #1 and #2. In this case, 30 different initial conditions for x i and y i are simulated and the maximum and minimum values of R are shown by the extremes in the filled area as a function of ε with α = 1.8. For ε 0.020, the networks depicts similar values of R . However, for 0.020 ε 0.053, this is changed, as we observe different values of R for the two cases. In this region the network has synchronization malleability: with only a difference between the sequence of inputs K i , the network can present R ≈ 0.88, for shuffling #1, or even R ≈ 0.03, for shuffling #2. For the region of high coupling strength, ε 0.053 the value of R is still different, but both networks depict phasesynchronized states. The results also show that the initial conditions of the systems do not seem to affect the dynamics of the network, since the filled area, representing variations due to different initializations, is only visible in a small region around ε ≈ 0.022. Therefore, synchronization malleability is observed for different inputs {K i }, but not for different initial conditions. We performed tests considering different values of σ and network size N and the synchronization malleability phenomenon can still be observed.
For a better visualization of the dynamical states exhibited by the network, Fig. 5 depicts raster plots (RP) of the spike times. Each black dot in the figure represents the time when a spike starts for the i-th neuron. We set α = 1.8 for all cases and consider different values of ε for shuffling #1 (first row) and shuffling #2 (second row). Regions indicated in the titles are the same ones defined in Fig. 3 (a) and Fig. 4 (d). Therefore, panels (a) and (f) are representative of non-synchronized states (ε = 0.005), where the spikes start at different times and there is no coherence in the neuronal activity. Panels (b) and (g) represent the case where ε = 0.021 and R is similar for both shuffles. In this case, one can observe diagonal structures and parts of the network with PS, which is in line with δR starting to depict high values (see Fig. 3

(b)).
A further increase in the coupling strength leads the networks to different dynamical states: panels (c) and (h) represent the case for ε = 0.031, where the network with shuffling #1 reaches PS ( R = 0.79) and horizontal structures can be observed, while the network with shuffling #2 does not reach PS ( R = 0.18) and, instead, diagonal structures can be observed. Instead of network PS, shuffling #2 has only groups of neurons with PS. For ε = 0.052, in panels (d) and (i), the dynamical difference between the networks increases: the network in (d) (shuffling #1) is phase-synchronized ( R = 0.88), but the one in (i) (shuffling #2) is not ( R = 0.03). In the latter case, one can notice the existence of locally horizontal structures in the raster plots, showing the existence of phase-synchronized communities that are not phase synchronized among themselves. This leads to the entire network having a small value of R . This situation is similar to the phenomenon of antiphase synchronization, where groups of oscillators depict, separately, synchronized characteristics, which cancel themselves globally [23,50]. At last, for the region of high coupling strength ε = 0.090, both networks depict phase-synchronized states and similar dynamical behaviors (panels (e) and (j)). The previous phenomena are observed for other shuffled input sequences and for other values of ε where the synchronization malleability is detected (see Fig. 3 (c)).
The results characterize the effect of the coupling strength variation for a fixed value of α. In this sense, as ε increases, one can observe a transition from non-phasesynchronized states to states where diagonals structures and local-phase synchronization are observed to, finally, phase-synchronized states. This kind of transition is observed for all networks with different shuffling processes for 1.5 α 2.0, which is the region in the parameter space where δR has high values (see Fig. 3 (c)). In order to analyze the spatiotemporal patterns, raster plots are plotted in the other panels of Fig. 6. In this case, the input values {K i } are given by shuffling #3 and #4. The values of α correspond to the regions 6, 7, and 8 defined in Fig. 3 (a) and Fig. 6 (a). The results shown in panels (b) and (c) (α = 1.0) indicate that both networks are in a phase-synchronized state, corroborating the high values of R observed in panel (a) for this region. However, for α = 1.8, the synchronization malleability phenomenon is again observed: panel (d) (shuffling #3) depicts a phase-synchronized state, with R = 0.92, while panel (e) (shuffling #4) depicts does not, with R = 0.05. At last, for higher values of α, the networks depict similar dynamical states without PS. Panels (f) and (g) (α = 2.5) show interesting spatiotemporal patterns with diagonal structures. This phenomenon can be understood as the effect of a stronger distance-dependence in the topology since higher values of α lead the effective coupling only with neighbors of each neuron. In this sense, coherence in the entire network is not obtained. Similar behavior is found in different models and coupling architecture when the coupling scheme has local characteristics [30,31,51]. To better understand why different shuffles of inputs {K i } lead to differences in the behavior of the network, we analyze its functional connectivity (FC). This is considered to be the matrix (MI) XY = M I(X, Y ) of mutual information between neurons, as described in section III C. Figure 7 shows the results of this procedure for the most noteworthy behaviors of the network.
The main importance of the FC analysis is the identification of the synchronization characteristics of the network since mutual information between pairs of neurons calculated in this way is proportional to their degree of PS [35]. In this context the results in Fig. 7 can be compared with the spatiotemporal pattern depicted by the raster plots of networks (Figs. 5 and 6): the horizontal lines observed in the FC are related with horizontal structures in the RP. It is clear then that, in cases where the network reaches PS, the horizontal structures in the FC are bigger. This means that neurons share more information and are more phase synchronized at both local (neighborhood) and global (long-range) levels.
The cases depicted in panels (a) and (b), with ε = 0.005, are simple: the network is desynchronized (small value of R ), which shows in the FC as mutual information being non-zero only on the main diagonal. A different scenario is observed in panels (c) and (d), where ε = 0.052. In panel (c), with shuffling #1, the FC is dominated by long horizontal and vertical structures (in yellow tones), indicating the sharing of information (PS) between both neighboring and distant neurons. Otherwise, for panel (d), with shuffling #2, horizontal and ver-tical structures are shorter, indicating a lesser degree of PS between neurons and, therefore, smaller R . In fact, the network in panel (c) has R = 0.88, while the one in (d) has R = 0.05. This is a case of synchronization malleability: some input sequences {K i } facilitate the formation of local and global structures, leading to PS, while others facilitate formation only of local structures, thus leading to some communities phase-synchronized within themselves, but not between themselves.
A further increase of the coupling strength to ε = 0.09 (panels (e) and (f)) leads both networks (with shuffling #1 and #2) to PS. The structures observed in panels (c) and (d) are intensified in this case, with values of the mutual information increasing. Both FCs now share information at local and global levels, meaning both networks are phase synchronized. However, panel (e) (shuffling #1) has longer horizontal and vertical structures than panel (f) (shuffling #2), indicating a higher degree of PS, which is indeed the case: R = 0.95 for (e) versus R = 0.82 for (f) (see Fig. 4 (d)).
To further illustrate the origin of the synchronization malleability, we show in panels (g) and (h) of Fig. 7 another example of different structures being formed due simply to a reordering of neurons' inputs (K i ). This is the extreme case depicted in Fig. 6 (d) and (e), where R = 0.92 and R = 0.05, respectively. Again, the case of global phase synchronization happens due to the formation of local and global structures, while the other happens due to the dominance of local structures.
At last, panels (i) and (j) are representative of states with diagonal structures and zig-zag states, observed for higher values of α (see panels (f) and (g) of Fig. 6). For both cases, higher values of mutual information between neurons are observed at diagonal lines in the FC, which correspond to the diagonal structures in the RP. In this case, the network is only frequency synchronized, but not phase synchronized. Furthermore, we note that, although R is similar to cases of very low coupling, like ε = 0.005 (for example, R ≈ 0.07 in panel (i) versus R ≈ 0.068 in (a)), the dynamics of the network is different, with a formation of various structures in this case.

V. CONCLUSIONS
Throughout this paper, we have analyzed a network composed of 525 spiking neurons simulated with the Chialvo map following a connection architecture described by a distance-dependent power-law scheme. We have shown that the interplay between coupling strength ε and power-law exponent α generates a great diversity of dynamical states. In networks with a weaker distance-dependence, a transition from desynchronized to phase synchronized states is observed with an increase in the coupling. By making the distance-dependence stronger, the network loses the phase-synchronized feature and there is a formation of new synchronization pat-terns, with locally phase-synchronized states and diagonal structures with only frequency synchronization. Similar behavior is found considering different local dynamics [46,48,52].
We have found a region in the parameter space α × ε where the networks depict a high level of sensitivity to changes in the neurons' inputs. In this case, very different dynamical states are observed for the same set of parameters (α, ε): networks can be either phase-synchronized or phase-desynchronized, depending on the ordering of the inputs, but not on the initial condition. By calculating the mutual information between the pairs of neurons in the network, we have obtained its functional connectivity and characterized its patterns of phase synchronization. We have seen that, for the malleable region, the sequences of inputs can either facilitate or hinder the formation of phase-synchronized structures. For sequences facilitating global structures, the network reaches phase synchronization. Otherwise, it does not, and phase synchronization is reached only between smaller groups of neurons.
Synchronization malleability was also found in a network of bursting neurons with coupling architecture following the Watts-Strogatz route, in which the network is taken from regular (local connections only) to smallworld to random by the rewiring of connections [49]. In this work, the phenomenon is observed when some local connections are changed by non-local ones, generating a coexistence of local (neighborhood) and global (longrange) topological effects. We suggest that a similar case also occurs in our results, with malleability occurring when the locality parameter α has intermediate values and leads to networks with a mix of local and global effectiveness.
Finally, this paper serves to characterize the phenomenon of synchronization malleability, for which there are still several open questions. Future works may study the mechanism behind this phenomenon, investigating the role of topology we proposed, and also look for it in other models and topologies.