Temporal-specific complexity of spiking patterns in spontaneous activity induced by a dual complex network structure

Temporal fluctuation of neural activity in the brain has an important function in optimal information processing. Spontaneous activity is a source of such fluctuation. The distribution of excitatory postsynaptic potentials (EPSPs) between cortical pyramidal neurons can follow a log-normal distribution. Recent studies have shown that networks connected by weak synapses exhibit characteristics of a random network, whereas networks connected by strong synapses have small-world characteristics of small path lengths and large cluster coefficients. To investigate the relationship between temporal complexity spontaneous activity and structural network duality in synaptic connections, we executed a simulation study using the leaky integrate-and-fire spiking neural network with log-normal synaptic weight distribution for the EPSPs and duality of synaptic connectivity, depending on synaptic weight. We conducted multiscale entropy analysis of the temporal spiking activity. Our simulation demonstrated that, when strong synaptic connections approach a small-world network, specific spiking patterns arise during irregular spatio-temporal spiking activity, and the complexity at the large temporal scale (i.e., slow frequency) is enhanced. Moreover, we confirmed through a surrogate data analysis that slow temporal dynamics reflect a deterministic process in the spiking neural networks. This modelling approach may improve the understanding of the spatio-temporal complex neural activity in the brain.

between excitatory/inhibitory neural populations [18][19][20] . Furthermore, the membrane potential behaviour under the threshold exhibits two particularly different states in spontaneous activity: the 'upstate' in which the membrane potential depolarizes near the threshold of spiking, and the 'downstate' in which the membrane potential hyperpolarizes and spikes cannot arise [18][19][20] .
Several studies over the last decade have used spiking neuron models to investigate the genesis of this spontaneous activity. Destexhe 21 showed that spontaneous activity may be produced by the existence of neurons with a low threshold. By using a different approach, investigators in another study showed that the structures of the neural network might induce spontaneous activity [22][23][24] . For example, Vogles & Abbott and Guo & Li showed that spontaneous activity is produced by the small-world property of synaptic connections characterized by small path lengths and large cluster coefficients and sparseness of random synaptic connections, respectively 22,23 . Of note, the model-which focuses on the distribution of the strength of cortical excitatory synaptic connections, as developed by Teramae et al. 24 -satisfies the physiological characteristics of spontaneous activity with high relevance.
In practical cerebral cortex, excitatory post-synaptic potentials (EPSPs), which is the increasing membrane potential by the pre-synaptic spikes in the excitatory synaptic connections, between cortical pyramidal neurons in a large majority of synapses exhibit potentials in the sub-millivolt range, whereas a sliver of synapses exhibit huge EPSPs (≳1.0 mV) 25,26 . In refs 25,26 , the distribution of EPSPs might follow a log-normal distribution. Regarding the mechanism that produce the spontaneous activity model proposed by Terame et al., spikes from most weak synapses (as noise) enhance spike transmission generated from a small number of strong synapses. An interpretation of this phenomenon could be that the mechanism of stochastic resonance induces spontaneous activity 24 . Furthermore, the EPSPs log-normal distribution may enhance the ability of associative memory recall 27 . It can also induce burst spiking, which has a vital function during hippocampal memory consolidation 28 .
Model-based studies of the spiking neural network with small-world properties have revealed that the complexity of neural activity is enhanced by the small-world network structure. Riecke et al. 29 demonstrated that complex spiking activity is induced by the small-world characteristics of excitable spiking neural networks through the rewiring process illustrated by the Watts and Strogatz model 30 . Shanahan 31 incorporated the physiological clustering features of the cortex into the spiking neural network by using the modular small-world network method 32 . Moreover, Watanabe et al. reported that the cerebral cortex exhibits a dual complex network structure, depending on the amount of EPSPs 33 . Weak synaptic networks and strong synaptic networks exhibit random network characteristics and small-world characteristics, respectively 33 . Under this circumstance, the functionalities for the duality of synaptic connectivity and the fluctuation of spontaneous activity have been of recent focus [34][35][36] .
In our previous study, the spiking neural network model for spontaneous activity 24 was expanded, and incorporated the previously described duality of synaptic connectivity, characterized by randomness and by small-worldness 37 . We attempted to demonstrate a tendency of enhancement of the complexity of spiking activity on a slow-temporal scale under the existing duality of synaptic connectivity 37 . However, physiological validation of the reproduced spontaneous activity remains to be achieved. In addition, the mechanism underlying this enhancement of complexity has not been elucidated.
Therefore, in this study, based on the outcome of our preliminary previous study 37 , we investigated the spontaneous activity in spiking neural networks using the duality of complex network structures, regarding these aforementioned points. First, we examined the physiological validity of the reproduced spontaneous activity at each duality level with respect to the low firing rate, coherence in spiking activity between excitatory/inhibitory neural populations, and coherence in the temporal spike series among neurons. Second, we analysed the multi-scale entropy (MSE) 38 of the time series of spontaneous activity and evaluated whether the obtained MSE reflected the dynamics of the spiking neural network using surrogate data analysis. Third, the mechanism for the enhancement of complexity was evaluated with a focus on the structures of spiking neural networks from the local network level and the entire topological level.

Methods
Spiking neural network with dual complex network structure. In this study, based on the spiking neural network with a log-normal EPSPs distribution proposed for producing the spontaneous activity 24 , the spiking neural network with a dual complex network structure, depending on synaptic weight, were constructed. In this network, we used the conductance-based leaky integrate-and-fire neuron model, to describe the membrane potential v(t): r thr ≥ → in which τ m is the decay constant of membrane. Also represented were the reversal potentials of the -amino-3-hydroxy-5-methyl-4-isoxazolepropionic acid (AMPA) receptor-mediated excitatory synaptic current, V E ; inhibitory synaptic current, V I ; leak current, V L ; and after-spike reset value of the v(t), V r . The I ex is the signal for external spikes triggering spontaneous activity, which is followed by A s δ(t − t ex ) mV. In this case, input time, t ex , is generated by the Poisson process, with a spiking rate of Λ Hz during [0:100] ms. The A s was set to V thr − V L + 1 to generate a spike from the resting state. The conductances for the excitatory synapse and inhibitory (represented by g E (t) [ms −1 ] and g I (t) [ms −1 ], respectively) are given by In this case, τ s is the synaptic decay constant. The s j is the spike time of the synaptic input from the j-th neuron. This spike has the synaptic delay, d j , and is enhanced by the excitatory synaptic weight, G E,j or diminished by the inhibitory synaptic weight, G I,j . In our study, parameter sets were used, as follows: V I = −80 mV; V L = −70 mV; V r = −60 mV; V thr = −50 mV; V E = 0 mV; τ m = 20 ms (excitatory neuron); τ m = 10 ms (inhibitory neuron); and τ s = 2 ms). In the numerical simulation, Eq. (3) was solved by the Euler method with a time-step size of Δt = 0.1 ms. The period of refractoriness was set to 1 ms. The synaptic delays in excitatory-to-excitatory connections were set to uniform random values [1:3] ms, and the synaptic delays for other connections were set to [0:2] ms. The sizes of this spiking neural network for excitatory neural population for and for inhibitory neurons were set to N E = 10000 and N I = 2000, respectively.
The EPSP amplitudes, V EPSP mV, which are increased membrane potentials resulting from the input spikes of excitatory synaptic connections, followed a log-normal distribution ( Fig. 1): 2 2 In this case, σ and the mode of the distribution were set to 1.0 and 0.2, respectively (μ − σ 2 = log0.2). The values V EPSP > 15 mV were declined as unrealistic values, and a new value was produced from the distribution.
Regarding the parameter to determine the topology for the excitatory-to-excitatory network, the number of synaptic connections was set to N . Synaptic connections with V EPSP ≤ 9 mV (i.e., weak synapses) were wired by using the Watts-Strogatz model 30 with a rewiring probability of 1.0, which corresponds a random network. However, synaptic connections with V EPSP > 9 mV (i.e., strong synapses) were wired with rewiring probability β from the ring type of topology. Each neuron in other synaptic connections was randomly connected; the coupling probability for excitatory connections and for inhibitory connections was 0.1 and 0.5, respectively.
The V EPSP is an observable value; therefore, these values must be converted into synaptic weight G E . In this conversion, we considered a simple case in which the spike input from a single excitatory synapse was applied to a post-synaptic neuron at t = 0 ms. The behaviour of the membrane potential, v(t), was as follows: We numerically solved Eqs (5) and (6). The relationship between V EPSP and G E can be clarified, as follows: The other synaptic weights were set to 0.018 (for excitatory-to-inhibitory connections), 0.002 (for inhibitory-to-excitatory connections), and 0.0025 (for inhibitory-to-inhibitory connections) 24 . The synaptic  24,26 . In this study, the spiking neural network was simulated by using Brian2 (https://brian2.readthedocs.io/en/2.0rc/index.html) 39 . www.nature.com/scientificreports www.nature.com/scientificreports/ Evaluation indexes. Index to measure neural activity. To measure spiking activity, the spiking rates in the excitatory neural population, r E Hz, and the inhibitory neural population, r I Hz, were used: In this instance, S E indicates the spike frequency in the bin in which the temporal width is 0.1 ms in excitatory neural populations, and S I is the spike frequency for inhibitory neural populations. A Gaussian-shaped window with temporal width 10 ms was applied to r E (t) and r I (t) to be smoothed in the time series. The period that has a high spiking rate r E > θ Hz, termed 'the activate period' in this study, was also used to quantify neural activity.
Index for the temporal complexity of neural activity. In this study, MSE 38 was utilised to quantify the complexity with time-scale dependency in the time series (r E ). A sample entropy for a stochastic variable {x 1 , x 2 , … x N } is given by the following equation: (i, j = 1, 2, …; i ≠ j). In this case, x i m is a vector with m dimension, given by = .
To evaluate temporal scale dependencies, a coarse-grained process for {x 1 , x 2 , …, x N } with the scale factor τ (τ = 1, 2, …) was applied in the MSE analysis: Against the coarse-grained time series, sample entropy h τ (r, m) was measured. For example, in the case of Owing to the dependency of h τ (r, m) on the scale factor τ, we evaluated the characteristics of the complexity in the time series of the excitatory firing rate r E (t) in 500 ≤ t ≤ 9500 [ms]. In our approach, we set m = 2 38 and set the width of scale at 1 ms. The value for r was set to 1.0 to fix the focusing patterns of r E , regardless of the value range. The MSE analysis was executed against 10 trials.
Indexes for topological features in synaptic connections. The topology of strong synaptic network (V EPSP > 9 mV), which consists of strong synapses as the edges and excitatory neurons as the nodes, were evaluated with regard to clustering, path length, and edge distribution. In this evaluation, a binarized adjacency matrix {a ij } (i, j = 1, 2, …, N E ) for strong synaptic connections was used.
The average clustering coefficient (CC) was defined by The average characteristic path length (PL) was defined as the average number of edges in the shortest path between two nodes in the network, as follows: Surrogate data analysis. We derived surrogate data using the iterative amplitude-adjusted Fourier transformed (IAAFT) surrogate data analysis 40 for the spiking rate, r E , to examine whether a nonlinear dynamic process was involved in the spiking. The iteration number was set to 50, and 10 IAAFT surrogate datasets were generated by different random seeds per original r E . These values of sample entropy h τ (r, m) were averaged and compared with the value from the original r E . To compare the sample entropies for the original data and those for IAAFT, we used the paired-sample t-test. The number of paired samples was set to 10. Two-tailed α levels of 10 −2 and 10 −3 were considered statistically significant.

Results
Characteristics of spontaneous activity and its physiological validation. The temporal characteristics of spiking activity was evaluated in the spiking neural network with a coexisting small-world network for strong synaptic connections and random network for weak synaptic connections under the conditions of different rewiring probabilities in strong synaptic connections: β = 1.0, 0.8, 0.6, 0.4, 0.2. The raster plot, time series of spiking rates of r E and r I , and power spectrum for spiking rates are shown in Fig. 2. In the β = 1.0, 0.8 case in which the topologies of a strong synaptic network (V EPSP > 9 mV) exhibits a nearly random network corresponding to the topologies of a weak synaptic network (V EPSP ≤ 9 mV), an irregular spatio-temporal spike pattern was shown in To examine whether the spiking activity at each β value was consistent with actual physiological spontaneous activity 18-20 , we evaluated the spiking rates, the correlation of temporal excitatory/inhibitory spiking activity, and the correlation among the temporal-spike series. Figures 3 and 4 depict the mean spiking rate, temporal correlation between excitatory/inhibitory neural populations, and the correlation among the temporal-spike series of excitatory neurons. Based on these results, we could observe that in all β cases, neurons exhibited low firing rates (r 10 E  Hz, r 100 I  Hz), high coherent spike transmission between specific neurons (the rate between highly coherent pairs with Pearson's correlation over 0.5. In addition, all pairs connected by the strong synapses were less than 4.2 × 10 −4 ), and high coherent spiking activity between excitatory/inhibitory neural populations (Pearson's correlation was approximately 0.98). This finding was consistent with the physiological findings [18][19][20] . These physiological conditions were primarily determined by lognormal-distribution of EPSPs 24 , whereas the dual complex network structure virtually did not affect them.
The MSE analysis for the time series of spontaneous activity. The complexity of temporal scale dependency for the time series of the spiking rate r E , corresponding to the results in Fig. 2, was analysed by MSE (given by Eqs 9 and 10). Figure 5 shows the profile of the sample entropy h τ of r E as function of temporal scale τ (corresponding frequency of 1/(10 −3 τ) Hz) for β = 1.0, 0.8, 0.6, 0.4, 0.2. In the β = 1.0, 0.8 case, dependency of h τ on temporal scale had a unimodal maximum peak at τ ≈ 20 (50 Hz). In smaller β cases (e.g., β = 0.6, 0.4, 0.2), the h τ increased at larger τ (i.e., small frequency) regions (τ ≳ 200 (5 Hz)). This characteristic was the same with the tendency at single trial shown in our previous preliminary results 37 .  www.nature.com/scientificreports www.nature.com/scientificreports/ The mechanism for the enhancement of complexity at a larger temporal scale was then evaluated. Figure 6 shows the activate period rate during the evaluation period as a function of β. The activate period rate decreased with increasing β value. Hence, the enhancement of complexity with the large temporal scale, as presented in Fig. 5, arose at a moderate activate period rate because of the intermittent appearance of the active state (see also Fig. 2).
Furthermore, we evaluated whether the temporal behaviour in spiking activity reflected a nonlinear dynamic process in the spiking neural networks. In the IAAFT process, the power spectrum was maintained, and the phase spectrum was randomised in the frequency space 40 . In comparison with the MSE profile for IAAFT surrogate data, a significant difference was observed at τ ≈ 200 (5 Hz) for cases of β = 0.6, 0.4, 0.2; the MSE profile under this condition reflected a nonlinear dynamic process (Fig. 7). This finding implied that the MSE profile at a large temporal scale reflected the characteristic of structure for the dual complex spiking network.
Relationship between spontaneous activity and network topology. We investigated the relationship between the MSE profile and the topology of a strong synaptic network (V EPSP > 9 mV) by clustering  www.nature.com/scientificreports www.nature.com/scientificreports/ coefficient and the mean shortest path length. Figure 8 shows the normalized clustering coefficient CC(β)/CC(0) (given by Eq. (12)) and mean shortest path length PL(β)/PL(0) (given by Eq. (13)), as a function of the rewiring probability β. Our results suggested that in the range 0.2 ≤ β ≤ 1.0 in which the β range for changing MSE profiles in Fig. 5), CC(β)/CC(0) decreased from 0.52 to 0. However, in this range of β, PL(β)/PL(0) was nearly zero 37 .
Our investigation confirmed that a neuron with high C i exhibits a tendency toward a high spiking rate in comparison with C i = 0 (see Fig. 9). Park et al. 41 reported that the local clustering properties increased the spiking rate in their small-world spiking neural network. Our results were congruent with these findings. With increasing β values, neurons with high C i decrease, leading to a decrease in the activate period for enhancing h τ at a large temporal scale.
In addition to the local clustering property, we investigated the distributions of the number of degrees D i around a node i and the spiking rate at each D i (see Fig. 10). The number of neurons with high D i increased with increasing β values, but the spiking rate and the activate period decreased (Figs 3 and 6). Therefore, neurons with a high degree of edge did not contribute to the generation of the activate state or the enhancement of h τ at a large temporal scale.

Discussion
In this study, we executed a simulation study using the leaky integrate-and-fire spiking neural network with log-normal synaptic weight distribution in EPSPs and with the duality of complex network consisting of a weak synaptic random network and a strong synaptic small-world network. We found that in cases in which the small-worldness of strong synaptic connections was enhanced, specific spiking patterns arose during temporal irregular spiking activity. The temporal scale profile of MSE for firing rates in the strong synaptic small-world network supported the enhancement of sample entropy at a larger temporal scale (slower frequency) than in cases in which strong synaptic connections approached those of random networks. These results were congruent with the single trial tendency shown in our previous preliminary results 37 . Furthermore, our previous work 37 mentioned that the clustering coefficient is a possible factor that produces enhancement at a large temporal scale. In this study, by the additional evaluation of local clustering characteristics, we confirmed that this enhancement of the sample entropy at a larger temporal scale was induced by the high spiking frequency of neurons with a high degree of clustering in strong synaptic connections. The IAAFT surrogate data analysis revealed that the sample entropy at the larger temporal scale reflected the dynamics the deterministic behaviour at the larger temporal scale might reflect the structural properties of spiking neural network, typified as the aforementioned strong synaptic clustering characteristic. The spiking activity, more importantly, also satisfied the physiological www.nature.com/scientificreports www.nature.com/scientificreports/ characteristics of spiking activity such as a low firing rate, high coherence in spiking activity between excitatory/ inhibitory neural populations, and high coherence in the temporal-spike series among a few neurons [18][19][20] . Based on these results, it can be interpreted that the enhancement of complexity at a large temporal scale with the characteristics of the physiological spontaneous activity can be induced by dynamics in the dual complex network structure as a physiological cortical network structure.
With regard to the emergence of various kinds of complex spatio-temporal spiking activity, initial studies 7-10 of spiking neural networks reported that random connections produced rich dynamics with highly irregular behaviour. Moreover, several studies applied physiological cortical structures to the random spiking network and showed that these expanded random network model can induce various spatio-temporal activities such as the intermittent transition between synchronous and asynchronous states [11][12][13][14] . Regarding the phenomena of emerging slow-temporal scale dynamics, Ostojic et al., Mastrogiuseppe & Ostojic and Wieland et al. reported that in a spiking neural network consisting of excitatory and inhibitory neural populations, irregular fluctuating activity, including slow-temporal scale dynamics, is induced in the strong strength regime of synaptic weight 12,13,15 . Moreover, Mart et al. have recently revealed that increasing symmetricity in random networks induces irregular dynamics with the slow-temporal scale 14 . They reported that, in a network with symmetricity (i.e., a portion of neurons have bidirectional connectivity), the neurons easily affect each other, and the neuron's activity has feedback. These interactions may induce the complex temporal behaviour with slow-temporal dynamics 14 . In our spiking neural network, in the regime of the strong clustering characteristic, the spikes of neurons belonging in same clusters easily affected each other and had feedback. Therefore, a similar effect in the random network with symmetricity 14 may arise in our spiking neural network.
Future directions and limitations of this study must be considered. During early development in the prenatal period, brain networks undergo enhancement of structural, ordered short-range connectivity. After birth, the short-range connectivity decreases, whereas long-range connectivity appears and strengthens, which leads to a decline in local clustering 42 . Furthermore, an increase in the complexity of the neural activity, as measured by electroencephalogram (EEG)/magnetoencephalography (MEG), is observed during development 43,44 . In particular, Hasegawa et al. 44 revealed that slower temporal complexity increases during early development. The changes that occur after birth can be interpreted as structural and functional integration. However, in Alzheimer's disease, the degradation of structural connectivity, including the progressive death of neurons, and the development of neurofibrillary tangles and senile plaques, induce a decline in the small-worldness property 45 . This structural degradation induces temporal-specific changes in the complexity of neural activity 4,46-52 . Our model focused on the duality of the cortical neural network structure, and did not consider physiological inter-regional connections. However, by considering region-specific subnetworks and inter-regional connections in this modelling approach may be a useful tool to reveal the relationships between the complexity of microscopic neural activity observed by neuroimaging modalities (e.g. EEG/MEG) and the topology of structural connectivity. Therefore, we plan to develop a spiking neural network consisting of region-specific subnetworks and inter-regional connections to evaluate the relationship between the complexity of neural activity and the topology of structural connectivity.
In conclusion, we conducted a simulation study by using a leaky integrate-and-fire spiking neural network that incorporated the duality of synaptic connectivity of a complex network that follows a log-normal synaptic weight distribution. It has been revealed that this duality has the function of generating complex neural activity, including complex dynamics on a large temporal scale. The future combination of this modelling approach with neuroimaging measures of the brain's networks and whole-brain network modelling may improve the understanding of the functionality of spatio-temporal complex neural activity.