Is chaos making a difference? Synchronization transitions in chaotic and nonchaotic neuronal networks

Chaotic dynamics of neural oscillations has been shown at the single neuron and network levels, both in experimental data and numerical simulations. Theoretical studies over the last twenty years have demonstrated an underlying role of chaos in neural systems. Nevertheless, whether chaotic neural oscillators make a significant contribution to relevant network behavior and whether the dynamical richness of neural networks are sensitive to the dynamics of isolated neurons, still remain open questions. We investigated transition dynamics of a medium-sized heterogeneous neural network of neurons connected by electrical coupling in a small world topology. We make use of an oscillatory neuron model (HB+Ih) that exhibits either chaotic or non-chaotic behavior at different combinations of conductance parameters. Measuring order parameter as a measure of synchrony, we find that the heterogeneity of firing rate and types of firing patterns make a greater contribution than chaos to the steepness of synchronization transition curve. We also show that chaotic dynamics of the isolated neurons do not always make a visible difference in process of network synchronization transitions. Moreover, the macroscopic chaos is observed regardless of the dynamics nature of the neurons. However, performing a Functional Connectivity Dynamics analysis, we show that chaotic nodes can promote what is known as the multi-stable behavior, where the network dynamically switches between a number of different semi-synchronized, metastable states.


Introduction
Over the past decades, a number of observations of chaos have been reported in the 2 analysis of time series from a variety of neural systems, ranging from the simplest to the 3 more complex [1,2]. It is generally accepted that the inherent instability of chaos in neurons [3,[13][14][15][16][17][18][19][20][21]. The first type of chaotic dynamics in neural systems is typically 14 accompanied by microscopic chaotic dynamics at the level of individual oscillators. The 15 presence of this chaos has been observed in networks of Hindmarsh-Rose neurons [8] and 16 biophysical conductance-based neurons [9][10][11][12]. The second type of chaotic firing pattern 17 is the synchronous chaos. Synchronous chaos has been demonstrated in networks of 18 both biophysical and non-biophysical neurons [3,13,15,17,[22][23][24], where neurons display 19 synchronous chaotic firing-rate fluctuations. The latter cases, the chaotic behavior is a 20 result of network connectivity, since isolated neurons do not display chaotic dynamics or 21 burst firing. More recently, a totally different mechanism showed that asynchronous 22 chaos, where neurons exhibit asynchronous chaotic firing-rate fluctuations, emerge 23 generically from balanced networks with multiple time scales synaptic dynamics [20]. 24 Different modeling approaches have been used to uncover important conditions for 25 observing these types of chaotic behavior (in particular, synchronous and asynchronous 26 chaos) in neural networks, such as the synaptic strength [25][26][27] in a network, 27 heterogeneity of the numbers of synapses and their synaptic strengths [28,29], and lately 28 the balance of excitation and inhibition [21] . The results obtained by Sompolinsky et 29 al. [25] showed that, when the synaptic strength is increased, neural networks displays a 30 highly heterogeneous chaotic state via a transition from an inactive state. Some other 31 studies demonstrated that chaotic behavior emerges in the presence of weak and strong 32 heterogeneities, for example a coupled heterogeneous population of neural oscillators 33 and the difference of their synaptic strengths [28][29][30]. Recently, Kadmon et al. [21] 34 highlighted the importance of balance of excitation and inhibition on a transition to 35 chaos in random neural networks. All these approaches identify the essential 36 mechanisms for generating chaos in neural networks. However, they give little insight 37 into whether chaotic neural oscillators make a significant contribution to relevant 38 network behavior, such as synchronization. In other words, whether the dynamical 39 richness of neural networks are sensitive to the dynamics of isolated neurons has not 40 been systematically studied yet. 41 To cope with this question, in the present paper we studied synchronization 42 transition in heterogeneous networks of interacting neurons. Here we make use of an 43 oscillatory neuron model (HB+I h ) that exhibits either chaotic or non-chaotic behavior 44 depending on biologically plausible parameter regions, as we showed in our previous 45 study [12]. We simulated small-world [31] neural networks consisting on a heterogeneous 46 population of HB+Ih neurons, connected by electrical synapses, and sampling their 47 parameters from either chaotic or non-chaotic regions of the parameter space. 48 Our first finding is that isolated chaotic neurons in networks do not always make a 49 visible difference in process of network synchronization transitions. The heterogeneity of 50 firing rate and the types of firing patterns make a greater contribution to the steepness 51 of the synchronization transition curve. Moreover, the macroscopic chaos is observed 52 regardless of the dynamics nature of the neurons. However, the results of Functional 53 Connectivity Dynamics (FCD) analysis show that chaotic nodes can promote what is 54 known as the multi-stable behavior, where the network dynamically switches between a 55 number of different semi-synchronized, metastable states. Finally, our results suggest 56 that chaotic dynamics of the isolated neurons is not a predictor of macroscopic chaos, 57 which can be a predictor of meta and multistability.

59
Single neuron dynamics 60 We use a parabolic bursting model inspired by the static firing patterns of cold 61 thermoreceptors, in which a slow sub-threshold oscillation is driven by a combination of 62 The membrane action potential of a HB+I h neuron follows the dynamics: where V is the membrane capacitance; I d , I r , I sd , I sr are depolarizing (Na V ), 70 repolarizing (K dr ), slow depolarizing (Na P / Ca T ) and slow repolarizing (K Ca ) currents, 71 respectively. I h stands for hyperpolarization-activated current, I l represents the leak 72 current, and lastly the term I syn is the synaptic current. Currents (except I syn ) are 73 defined as: where a i is an activation term that represents the open probability of the channels 75 (a l ≡ 1), with the exception of a sr that represents intracellular Calcium concentration.

76
Parameter g i is the maximal conductance density, E i is the reversal potential and the 77 function ρ(T ) is a temperature dependent scale factor for the current.

78
The activation terms a r , a sd and a h follow the differential equations: where On the other hand, a sr follows Finally, The function φ(T ) is a temperature factor for channel kinetics. The In the simulation, we vary the maximal conductance density g sd , g sr and g h values. 86 Unless stated otherwise, the parameters used are given in Table 1.

88
The synaptic input current into neuron k is given by: In this article, the current I kl between neuron k and l is modeled as a gap junction 90 (electrical synapse): Where g kl is the conductance (coupling strength) of synapse from cell k to cell l. In 92 this work we selected a uniform value g kl = g for all connections within the neural 93 networks, and simulations were performed with different values of g in order to observe 94 the transition to total network synchrony.

95
Structural connectivity matrix. 96 We define the connectivity matrix by C kl = 1 if the neuron k is connected to neuron l 97 and C kl = 0 otherwise. We employed a Newman-Watts small world topology [31], 98 implemented as two basic steps of the standard algorithm [33][34][35]: (1). create a ring 99 lattice with N nodes of mean degree 2K. Each node is connected to its K nearest 100 neighbors on either side. (2). For each node in the graph, add an extra edge with 101 probability p, to a randomly selected node. The added edge cannot be a duplicate or  The method for establishing whether a system is chaotic or not is to use the Lyapunov 108 exponents. In particular, Maximal Lyapunov exponent (MLE) greater than zero, hence, 109 is widely used as an indicator of chaos [36][37][38][39]. We calculated MLEs from trajectories in 110 the full variable space, following a standard numerical method based on that of Sprott 111 (2003) [36] (see also Jones et al., 2009) [37].

113
The voltage trajectory of each neuron was low-pass filtered (50Hz) and a continuous 114 Wavelet transform [40][41][42] was applied to determine in one step the predominant 115 frequency and the instantaneous phase at that frequency. We use the complex Morlet 116 wavelet as mother wavelet function to calculate instantaneous phase. 117 We describe global dynamical behavior of the neural networks using the mean and 118 the standard deviation of the order parameter amplitude over a time-course, which 119 indicate respectively the global synchrony and the global metastability of the 120 system [43,44]. The order parameter [45,46], R, describes the global level of phase 121 synchrony in a system of N oscillators, given by: Where ϕ k (t) is the phase of oscillator k at time t, f N = φ c (t) denotes the average 123 of f over all k in networks, | | is absolute value and f t is the average in time.
corresponds to the maximally asynchronous (disordered) state, whereas R = 1 125 represents the state where all oscillators are completely synchronized (phase synchrony 126 state). The global metastability χ of neural networks is given by: Metastability is zero if the system is either completely synchronized or completely 128 desynchronized during the full simulation-a high value is present only when periods of 129 coherence alternate with periods of incoherence [43], where ∆t of Eq. (12) is time 130 windows to quantify the global metastability.

132
A functional connectivity (FC) matrix was calculated by the pair-wise phase synchrony 133 (order parameter). This was done in a series of M overlapping time windows Where t corresponds to all times inside window m. The values in the lower triangle 136 of FC, discarding the diagonal and the values adjacent to it, were vectorized and a 137 correlation matrix is calculated between the vectors. This constitutes the Functional 138 Connectivity Dynamics (FCD) matrix [47,48]: Where cov (X, Y ) is the covariance between vectors X and Y , and σ(X) is the

Synchronization transitions with parameters drawn from fixed-size regions 152
Our main goal is to study how the dynamics of isolated neural oscillators can propagate 153 to the network level, in terms of relevant behaviors. To do this, we use a model of 154 neural oscillator that can display either chaotic or non-chaotic behavior depending on 155 the parameters (Fig 1, Fig 2A, Xu et al 2017 [12]). We simulated networks of 50 156 neurons connected by electrical coupling (gap junctions) in a small world topology.

157
When drawing the g sd and g h parameters, we selected different regions of the parameter 158 space that made them behave as either chaotic or non-chaotic oscillators, while 159 maintaining a similar average firing rate (Fig 2A). Then, the inter-cellular conductance 160 g was varied from 0 to 1 in order to evidence the transition from asynchrony to 161 complete synchronization. oscillators indeed impacts the network behavior. In networks of non-chaotic oscillators, 169 when g is between 10 −4 and 10 −2 , we observe that the networks becomes chaotic while 170 transition from asynchrony (low R) to synchrony (high R). However, when phase 171 synchronization is reached, chaos is lost. Thus, only asynchronous chaos is observed. In 172 the case of chaotic neurons, we find that the networks always exhibits chaotic behavior 173 at a wide range (with g from 0 to 1) of synaptic coupling. In that case, asynchronous 174 and synchronous chaos both are observed. When the same analysis was applied to the 175 second pair of parameter regions, we observe a strange behavior of the non-chaotic 176 oscillators, with a non-smooth transition associated to a higher metastability and 177 network MLE. An inspection of the firing patterns (Fig 2A, right) reveals that the 178 regions labeled as '2' contains transitions between different firing patterns: tonic regular 179 to bursting for non-chaotic, and skipping to bursting for the chaotic region. This made 180 us think that the 'kink' observed in the curves of non-chaotic oscillators, was due to a 181 transition between firing patterns occurring in the network. It can be seen in Fig 3A   182 that as synaptic coupling is increased in networks of non-chaotic neurons, bursting firing 183 pattern dissapear, while in chaotic neurons they are increased. Moreover, in non-chaotic 184 neurons there is a rebound of bursting firing patterns at g = 0.0433. This finding 185 suggests that the dramatical changes in firing patterns induce the non-smooth 186 transition shown in Fig 2B,(bottom). non-chaotic. Thus, the steeper transition to synchrony and lower metastability observed 192 with non-chaotic oscillators could be due to a more homogeneous nature of the network, 193 rather than the chaotic oscillation by itself. Although this is already an effect of chaos, 194 it is of dubious biological relevance because neurons will never control their levels of 195 channel expression within a fixed range, as we did here. If any, neurons control for 196 function, and a simple approximation to this is to consider that they try to maintain a 197 certain average firing rate with whatever ion channel density relationship that can 198 attain it. Thus, we developed a parameter sampling procedure that replicated, for both 199 chaotic and non-chaotic populations, a similar firing rate distribution rather than the 200 mean. We also shifted to another parameter sub-space (g sd /g sr ) to take advantage of 201 the complete absence of chaos when g h = 0 (Xu et al [12]). Nevertheless, the 202 simulations that follow were also performed in the g sd /g h parameter subspace with very 203 similar results (see Supplemental Material S1)  [12] and Fig 4B). The desired region of this case in 215 firing rate is shown in the Fig 4B (bottom) . The histogram of firing rates in each 216 selected set of parameters (Fig 4C), shows that they have the same distribution of firing 217 rate. The same operation was used to select parameter sets with higher firing rate (from 218 7.0 to 9.5 Hz, not shown).  236 Finally, we measured the ability of our network models to display multi-stable behavior 237 by characterizing their functional connectivity dynamics (FCD). This analysis is being 238 extensively applied to fMRI and M/EEG recordings [47][48][49] and is explained in Fig 6   239 and Methods. Briefly, the series is divided in overlapping time windows and for each   Fig 7A below each FCD) are also useful in detecting the three 254 situations described. As a rough measure of multi-stability, we took the variance of the 255 FCD values (outside the diagonal) and plotted them against the synaptic conductance, 256 averaging several simulations with different seed for the random connectivity matrix 257 (Fig 7B). In the 3.0 to 4.5 firing rate range (for the isolated oscillators), it is clear that 258 chaotic nodes produce networks with higher multi-stability than both non-chaotic and 259 NoIh nodes. Moreover, the g range in which the multi-stable behavior is observed is 260 wider. In the 7.0 to 9.5 firing rate range, the variance of the FCD is not higher for 261 chaotic nodes, however the g range for multi-stability is still wider. This shows that

266
In this work, we investigated how a complex node dynamics can affect the 267 synchronization behavior of a medium-sized heterogeneous neural network. We maximum MLEs were found. However this must be taken with caution as the MLE is 288 not necessarily a quantitative measure of chaos [51].

289
The chaotic nature of the isolated oscillators did not always convey to network chaos 290 in a direct or predictable fashion. More specifically, our networks always showed chaotic 291 behavior at some g values regardless of the dynamics of the isolated nodes. This is not 292 surprising, as chaotic behavior arises in networks of very simple units and seems to 293 depend more strongly on other factors such as synaptic weights and network 294 topology [25][26][27]52]. Moreover, just the high-dimensionality of the systems seems to be 295 enough to assure that chaos will emerge under some conditions, for example, the 296 quasiperiodic route to chaos in high-dimensional systems [53].

297
Assessing chaos in a large network or in a high-dimension system can be a difficult 298 task. It is now generally accepted that a unique intrinsic and observable signature of 299 systems exhibiting deterministic chaos is a fluctuating power spectrum with an 300 exponential frequency dependency [54]. Thus, some studies introduced the broad power 301 spectrum to characterize the chaos of networks [13,55]. Here we use the most popular 302 and directly method of maximal Lyapunov exponent (MLE) to quantify chaos on the 303 level of networks in the way as we did for single cells [12,36,37]. As usual, we define the 304 state of the network as chaotic if the MLE is greater than zero.

305
Macroscopic chaos mentioned above arises from the network's global properties, the 306 propensity of isolated neurons to oscillate, the nature of synaptic dynamics, or a 307 mixture of the them, as shown in earlier works [25][26][27][28][29]. In this paper, the focus is 308 different. First, finding both asynchronous and synchronous chaos in the same network, 309 only by changing the synaptic strength, is new. Secondly, the route from asynchronous 310 to synchronous chaos in networks of chaotic and non-chaotic oscillations has a slight 311 difference and has not been found in previous studies. Specifically, network switches 312 directly from asynchronous to synchronous chaos in networks of chaotic neurons, while 313 networks of non-chaotic neurons usually can go through four phases of network state, 314 that are asynchronous activity, asynchronous chaos, then again asynchronous activity 315 and lastly synchronous chaos.

316
While discussing chaos in neural systems we have used completely deterministic 317 dynamics. Random variables were used to define network connectivity and node 318 parameters, but the time evolution of the networks and nodes was calculated in the 319 absence of noise. However, neural systems are subject to a number of noise sources, 320 being the most important the stochastic opening an closing of ion channels and synaptic 321 variability [56]. How the synchronization transitions and meta/multi-stable behavior 322 will emerge in a noisy system remains to be studied, and it will be interesting to assess 323 how much the dynamics introduced by chaos can prevail in the presence of noise.

324
In summary, we have shown that chaotic neural oscillators can make a significant   what has been plotted of main text (shown in Fig 5).