Introduction

Recent progress in studies using neuroimaging modalities has elucidated the structure of the whole network of the brain, which is composed of individual regions and connections, called a connectome1,2. In communications among neuronal populations or among brain regions, feedback loops at multiple hierarchical levels of cortical processing can produce spatio-temporal complex behaviours that correspond to brain activity influenced by recurrent neural processes3,4. In a recent modelling approach, which utilises network structures such as topological features and coupling strength estimated by neuroimaging modalities, whole-brain spiking neural networks have been constructed and reproduced complex neural behaviours5,6. While, despite the intra-regional cortical networks level without these hierarchical topological structures, the neural network exhibits various types of complex spatio-temporal behaviours7,8,9,10,11,12,13,14,15. To describe these behaviours, randomly connected network models, which are approximations for complex cortical networks, were initially proposed, and produced highly irregular behaviours through mutual interactions7,8,9,10. Moreover, these models, which include physiological cortical network structures such as electrical coupling, excitatory and inhibitory neural populations, spiking neurons, and symmetricity in random networks, can reproduce various types of spatio-temporal neural activity with higher physiological validations in comparison with random networks11,12,13,14,15.

With regard to the typical behaviour of the neural activity in the cortex, this network represents the sustaining electrical dynamic fluctuation, despite the nonexistence of sensory stimulation, which is called spontaneous activity16,17. This neural activity has several major characteristics: irregular low-frequency neural spiking (≈1 Hz) with high coherent spike transmission between specific pairs of neurons, and high coherent spiking activity between excitatory/inhibitory neural populations18,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 arise18,19,20.

Several studies over the last decade have used spiking neuron models to investigate the genesis of this spontaneous activity. Destexhe21 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 activity22,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, respectively22,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 refs25,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 activity24. Furthermore, the EPSPs log-normal distribution may enhance the ability of associative memory recall27. It can also induce burst spiking, which has a vital function during hippocampal memory consolidation28.

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 model30. Shanahan31 incorporated the physiological clustering features of the cortex into the spiking neural network by using the modular small-world network method32. Moreover, Watanabe et al. reported that the cerebral cortex exhibits a dual complex network structure, depending on the amount of EPSPs33. Weak synaptic networks and strong synaptic networks exhibit random network characteristics and small-world characteristics, respectively33. Under this circumstance, the functionalities for the duality of synaptic connectivity and the fluctuation of spontaneous activity have been of recent focus34,35,36.

In our previous study, the spiking neural network model for spontaneous activity24 was expanded, and incorporated the previously described duality of synaptic connectivity, characterized by randomness and by small-worldness37. 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 connectivity37. 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 study37, 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 activity24, 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):

$$\frac{dv}{dt}=-\frac{1}{{\tau }_{m}}(v-{V}_{L})-{g}_{E}(v-{V}_{E})-{g}_{I}(v-{V}_{I})+{I}_{{\rm{ex}}},$$
(1)
$${\rm{if}}\,v\ge {V}_{{\rm{thr}}}\,{\rm{mV}},\,{\rm{then}}\,v(t)\to {V}_{r},$$
(2)

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, VE; inhibitory synaptic current, VI; leak current, VL; and after-spike reset value of the v(t), Vr. The Iex is the signal for external spikes triggering spontaneous activity, which is followed by Asδ(t − tex) mV. In this case, input time, tex, is generated by the Poisson process, with a spiking rate of Λ Hz during [0:100] ms. The As was set to Vthr − VL + 1 to generate a spike from the resting state. The conductances for the excitatory synapse and inhibitory (represented by gE(t) [ms−1] and gI(t) [ms−1], respectively) are given by

$$\frac{d{g}_{X}}{dt}=-\frac{{g}_{X}}{{\tau }_{s}}+\sum _{j}\,{G}_{X,j}\sum _{{s}_{j}}\,\delta (t-{s}_{j}-{d}_{j}),\,X=E,I.$$
(3)

In this case, τs is the synaptic decay constant. The sj is the spike time of the synaptic input from the j-th neuron. This spike has the synaptic delay, dj, and is enhanced by the excitatory synaptic weight, GE,j or diminished by the inhibitory synaptic weight, GI,j. In our study, parameter sets were used, as follows: VI = −80 mV; VL = −70 mV; Vr = −60 mV; Vthr = −50 mV; VE = 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 NE = 10000 and NI = 2000, respectively.

The EPSP amplitudes, VEPSP mV, which are increased membrane potentials resulting from the input spikes of excitatory synaptic connections, followed a log-normal distribution (Fig. 1):

$$p(x)=\frac{\exp [-\,{(\log x-\mu )}^{2}/2{\sigma }^{2}]}{\sqrt{2\pi }\sigma x}.$$
(4)
Figure 1
figure 1

Synaptic weight distribution among excitatory neurons (upper: single logarithmic plot; lower: double logarithmic plot). The dotted line represents the threshold dividing the network structures. EPSP: excitatory postsynaptic potential. (σ = 1.0, μ − σ2 = log0.2).

Full size image

In this case, σ and the mode of the distribution were set to 1.0 and 0.2, respectively (μ − σ2 = log0.2). The values VEPSP > 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 \(0.1{N}_{E}^{2}\). Synaptic connections with VEPSP ≤ 9 mV (i.e., weak synapses) were wired by using the Watts-Strogatz model30 with a rewiring probability of 1.0, which corresponds a random network. However, synaptic connections with VEPSP > 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 VEPSP is an observable value; therefore, these values must be converted into synaptic weight GE. 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:

$$\frac{dv(t)}{dt}=-\frac{1}{{\tau }_{m}}(v(t)-{V}_{L})-{g}_{E}(t)(v(t)-{V}_{E}),$$
(5)
$$\frac{d{g}_{E}(t)}{dt}=-\frac{{g}_{E}(t)}{{\tau }_{s}}+{G}_{E}\delta (t).$$
(6)

We numerically solved Eqs (5) and (6). The relationship between VEPSP and GE can be clarified, as follows:

$${G}_{E}={V}_{{\rm{EPSP}}}/100.$$
(7)

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 transmissions in the excitatory-to-excitatory synaptic connections failed, based on the failing rate, depending on the amplitude of the EPSP: \({P}_{E}=\frac{a}{a+{V}_{{\rm{EPSP}}}}\) (a = 0.1 mV)24,26. In this study, the spiking neural network was simulated by using Brian2 (https://brian2.readthedocs.io/en/2.0rc/index.html)39.

Evaluation indexes

Index to measure neural activity

To measure spiking activity, the spiking rates in the excitatory neural population, rE Hz, and the inhibitory neural population, rI Hz, were used:

$${r}_{X}(t)={10}^{3}\frac{{S}_{X}(t)}{\Delta t{N}_{X}}\,X=E,I.$$
(8)

In this instance, SE indicates the spike frequency in the bin in which the temporal width is 0.1 ms in excitatory neural populations, and SI is the spike frequency for inhibitory neural populations. A Gaussian-shaped window with temporal width 10 ms was applied to rE(t) and rI(t) to be smoothed in the time series. The period that has a high spiking rate rE > θ 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, MSE38 was utilised to quantify the complexity with time-scale dependency in the time series (rE). A sample entropy for a stochastic variable {x1, x2, … xN} is given by the following equation:

$$h(r,m)=-\,\log \,\frac{{C}_{m+1}(r)}{{C}_{m}(r)},$$
(9)

in which Cm(r) represents the probability of satisfying the condition with \(|{{\bf{x}}}_{i}^{m}-{{\bf{x}}}_{j}^{m}| < r\) (i, j = 1, 2, …; i ≠ j). In this case, \({{\bf{x}}}_{i}^{m}\) is a vector with m dimension, given by

$${{\bf{x}}}_{i}^{m}=\{{x}_{i},{x}_{i+1},\cdots ,{x}_{i+m}\}.$$
(10)

To evaluate temporal scale dependencies, a coarse-grained process for {x1, x2, …, xN} with the scale factor τ (τ = 1, 2, …) was applied in the MSE analysis:

$${y}_{j}^{(\tau )}=\frac{1}{\tau }\mathop{\sum }\limits_{i=(j-1)\tau +1}^{{j}_{\tau }}\,{x}_{i}.\,(1\le j\le N/\tau )$$
(11)

Against the coarse-grained time series, sample entropy hτ(r, m) was measured. For example, in the case of τ = 2, \(\{{y}_{j}^{(2)}\}=\{\frac{{x}_{1}+{x}_{2}}{2},\frac{{x}_{3}+{x}_{4}}{2},\cdots ,\frac{{x}_{N-1}+{x}_{N}}{2}\}\). 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 rE(t) in 500 ≤ t ≤ 9500 [ms]. In our approach, we set m = 238 and set the width of scale at 1 ms. The value for r was set to 1.0 to fix the focusing patterns of rE, 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 (VEPSP > 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 {aij} (i, j = 1, 2, …, NE) for strong synaptic connections was used.

The average clustering coefficient (CC) was defined by

$$CC=\frac{1}{{N}_{E}}\sum _{i\in {N}_{E}}\,\frac{2{C}_{i}}{{k}_{i}({k}_{i}-1)},$$
(12)

in which ki and Ci represent the degree of a node i, given by \({k}_{i}\sum _{j\in {N}_{E}}\,{a}_{ij}\), and the number of triangles around node i, respectively.

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:

$$PL=\frac{1}{{N}_{E}({N}_{E}-1)}\sum _{i,j\in {N}_{E},j\ne i}\,{d}_{ij},$$
(13)

in which the matrix for the shortest path length between nodes i and j are represented by {dij}. In the Watts–Strogatz model, the distribution of the degree of edges changes, depending on β. Therefore, we assessed the degree of node i using Di. To observe the dependency of topological changes on rewiring properties β in the Watts–Strogatz model, the normalized average clustering coefficient was used, based on the β = 0 case [CC(β)/CC(0) and PL(β)/PL(0)].

Surrogate data analysis

We derived surrogate data using the iterative amplitude-adjusted Fourier transformed (IAAFT) surrogate data analysis40 for the spiking rate, rE, 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 rE. These values of sample entropy hτ(r, m) were averaged and compared with the value from the original rE. 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 rE and rI, 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 (VEPSP > 9 mV) exhibits a nearly random network corresponding to the topologies of a weak synaptic network (VEPSP ≤ 9 mV), an irregular spatio-temporal spike pattern was shown in the raster plot. The time series of spiking rates exhibited irregular oscillations around rE ≈ 4.0 Hz and rI ≈ 40 Hz. Moreover, constant spike patterns arose in a specific duration in cases of a smaller β value (β = 0.6, 0.4, 0.2) such as [6500:6600] ms for β = 0.6, [6600:7250] ms for β = 0.4, and [6000:8000] for β = 0.2 (see the raster plot). The spiking rate was increased for rE, rI, compared with cases in which β = 1.0, 0.8 (rE ≈ 5.0 Hz and rI ≈ 70 Hz for β = 0.6; rE ≈ 6.0 Hz and rI ≈ 100 Hz for β = 0.4; and rE ≈ 7.0 Hz and rI ≈ 110 Hz for β = 0.2). In addition, the power spectrum of the spiking rates for β = 1.0, 0.8, 0.6, 0.4, 0.2 was mostly distributed in the region \(\lesssim 40\) Hz (power/Freq ≳ −40 [dB/Hz]).

Figure 2
figure 2

Rewiring probability β dependency of spiking activity in a spiking neural network with small-world network for strong synaptic connections and random network for weak synaptic connections (β = 1.0, 0.8, 0.6, 0.4, 0.2). The left, middle, and right images represent the raster plot, spiking rate time-series for rE (upper) and rI (lower), and their power spectrum, respectively. In the β = 1.0, 0.8 case, an irregular spatio-temporal spike pattern was observed. With decreasing β values (β = 0.6, 0.4, 0.2), constant spike patterns with a higher spiking rate [rE ≈ 5.0 Hz, (rI ≈ 70 Hz) in the β = 0.6 case, rE ≈ 6.0 Hz (rI ≈ 100 Hz) in the β = 0.4 case, and rE ≈ 7.0 Hz (rI ≈ 110 Hz) in the β = 0.2 case] were maintained in a specific duration.

Full size image

To examine whether the spiking activity at each β value was consistent with actual physiological spontaneous activity18,19,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}_{E}\lesssim 10\) Hz, \({r}_{I}\lesssim 100\) 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 findings18,19,20. These physiological conditions were primarily determined by lognormal-distribution of EPSPs24, whereas the dual complex network structure virtually did not affect them.

Figure 3
figure 3

The mean spiking rate in excitatory/inhibitory neural populations (upper) and the correlation between the time series of the excitatory spiking rate rE and that of the inhibitory spiking rate rI (lower). For all β cases, neurons exhibit low firing rates (\({r}_{E}\lesssim 0\) Hz, \({r}_{I}\lesssim 100\) Hz) and high coherent spiking activity between the excitatory and inhibitory neural populations.

Full size image
Figure 4
figure 4

Neuron pairs satisfying the high correlation among the temporal spike series. The black dots indicate neuron pairs with Pearson’s correlation over 0.5. High coherent spike transmission between specific neurons. Rates among pairs with Pearson’s correlation over 0.5 in all pairs connected by strong synapses are 2.6 × 10−7 (β = 1.0), 2.8 × 10−7 (β = 0.8), 4.9 × 10−5 (β = 0.6), 4.2 × 10−4 (β = 0.4), and 2.4 × 10−4 (β = 0.2).

Full size image

The MSE analysis for the time series of spontaneous activity

The complexity of temporal scale dependency for the time series of the spiking rate rE, 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 rE 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 results37.

Figure 5
figure 5

Profile of multiscale entropy as a function of temporal scale τ for the time series of the excitatory spiking rate rE (the width of the temporal scale is 1 ms). The mean and standard deviation between trials are represented by the solid and shaded areas, respectively. In the β = 1.0, 0.8 case, hτ exhibits a peak at τ ≈ 20 (50 Hz). With decreasing β values such as β = 0.6, 0.4, 0.2, the hτ in the larger τ (i.e., small frequency) regions increases (τ ≳ 200 (5 Hz)).

Full size image

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).

Figure 6
figure 6

Dependence of the average of the activate period for the excitatory neural population (i.e., satisfies rE > θ Hz) on rewiring probability β. The activate period rate decreases with increasing β value.

Full size image

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 space40. 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.

Figure 7
figure 7

Comparison between multiscale entropy profiles of the original rE and those of the IAAFT surrogate data. The cases in which differences in hτ are statically significant are represented by blue (+p < 10−2) and red (*p < 10−3). A significant difference between the original/IAAFT ones is observed at τ ≈ 200 (5 Hz) in cases of β = 0.6, 0.4, 0.2; the MSE profile at large temporal scales reflects a nonlinear dynamic process. Org: original; IAFFT: Iterative amplitude-adjusted Fourier transformed.

Full size image

Relationship between spontaneous activity and network topology

We investigated the relationship between the MSE profile and the topology of a strong synaptic network (VEPSP > 9 mV) by clustering 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 zero37.

Figure 8
figure 8

In the strong synaptic network, the normalised clustering coefficient CC(β)/CC(0) and the mean shortest path length PL(β)/PL(0), as a function of the rewiring probability β. In the range 0.2 ≤ β ≤ 1.0 (which corresponds to the β range for changing MSE profiles in Fig. 5), CC(β)/CC(0) decreases from 0.52 to 0, whereas PL(β)/PL(0) remains at approximately 0.

Full size image

Our investigation confirmed that a neuron with high Ci exhibits a tendency toward a high spiking rate in comparison with Ci = 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 Ci decrease, leading to a decrease in the activate period for enhancing hτ at a large temporal scale.

Figure 9
figure 9

Relationship between spiking activity and clustering structures. Probability distribution of the number of triangles around node i: Ci (upper) and dependence of the temporal average spiking rate of the excitatory neural population on Ci (lower). A neuron with high Ci exhibits a tendency toward a high spiking rate, compared with the case of Ci = 0. Ci: number of triangles around node i. Ave.: average; Prob.: probability; Dist.: distribution.

Full size image

In addition to the local clustering property, we investigated the distributions of the number of degrees Di around a node i and the spiking rate at each Di (see Fig. 10). The number of neurons with high Di 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.

Figure 10
figure 10

Probability distribution of the number of degrees of node i: Di. The number of neurons with high Di increased with increasing β values, but the spiking rate and the activate period decreased. Prob.: probability; Dist.: distribution.

Full size image

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 results37. Furthermore, our previous work37 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 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 neurons18,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 studies7,8,9,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 states11,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 weight12,13,15. Moreover, Mart et al. have recently revealed that increasing symmetricity in random networks induces irregular dynamics with the slow-temporal scale14. 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 dynamics14. 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 symmetricity14 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 clustering42. Furthermore, an increase in the complexity of the neural activity, as measured by electroencephalogram (EEG)/magnetoencephalography (MEG), is observed during development43,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 property45. This structural degradation induces temporal-specific changes in the complexity of neural activity4,46,47,48,49,50,51,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.