Overrepresented connectivity patterns.

To demonstrate the presence of overrepresented pair and triad motifs in our model as found in the layer 5 of rat visual cortex [5], we first compare the number of bidirectional connections against Erdős-Rényi (E-R) random networks. We find significantly more bidirectional connections than expected in E-R networks (mean ± SEM: 1.94 ± 0.02, p < 0.001, one sample t-test), indicating that the bidirectional connections are overrepresented. A triad motif can be arranged in thirteen ways with all three neurons connected to each other (Fig 2A). As in [5], we compare the numbers of these motifs in our network against those in a null model, which is a random network with the unidirectional and bidirectional connection probabilities specified separately to match the original network. As shown in Fig 2A, the triad motifs in our model exhibit two salient features. First, for the most overrepresented motif, No.13 (the fully connected motif), its count relative to the null model (4.96 ± 1.86; p < 0.001, one sample t-test) is dominantly higher than the other motifs; this feature of the fully connected motif in our network model is consistent with measured data [5]. Second, the majority of the motifs that form a closed loop (ignoring direction) are overrepresented (except motif No.7 and No.11, Fig 2A); this is also consistent with experimental observations [5]. The existence of such closed loops in a triad motif can be characterized by its directed clustering coefficient (CC). The CC of a node (neuron) indicates the tendency of the node’s neighbors to cluster together [27], with its value ranging from 0 (none of the neighbors are connected to each other) to 1 (all neighbors are mutually reciprocally connected); the CC of a motif is obtained by averaging the CC of its three nodes. As shown in Fig 2A, the motifs with a closed loop have non-zero CCs and vice versa. As will be shown in the later sections, the CC is a structural property that has significant impact on neural dynamics.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Emergent higher-order connectivity features.

(A) Triad motif counts. Within each trial, the motif counts are averaged from 1000 randomly sampled quadruplets using the brain connectivity toolbox [28] and 100 control networks are generated to normalize the motif counts. The error bars show one SEM from 10 trials. The mean clustering coefficient is shown above each motif, with the non-zero ones highlighted in red. The asterisks denote the relative motif counts that are statistically significantly different from one (red asterisks p < 0.001, black asterisks p < 0.034; Holm step-down adjusted p-values [29]). (B) Left y-axis (blue): normalized rich-club coefficients Φ(k), which measures the extent to which the neurons with total degrees > k in a network are over-connected to each other. Within each trial, 100 control networks are generated to normalize the coefficients. The coefficient peaks at 1.60, which is significantly higher than 1 as given by the control networks (p < 0.001, one-sample t-test). Right y-axis (red): dynamical importance (averaged across neurons with the same total degree k within each trial), defined as the fractional change (%) in the largest eigenvalue of the adjacency matrix of the network upon removal of the neuron from it. The solid lines show the averages over 10 trials and the shaded areas show one SEM. Data for k < 200 is omitted from the plot because neurons with total degree k < 200 are sparse. (C) The rich-club curves with different common neighbor coefficients a_Γ. Inset: the maximum rich-club coefficient between k = 500 and k = 600 are positively correlated with the common neighbor coefficient. Both motif and rich-club results are calculated for a circular region in the network containing 2000 excitatory neurons to avoid artefacts from the periodic boundaries.

https://doi.org/10.1371/journal.pcbi.1006902.g002

In our model, the overrepresentation of the triad motifs with non-zero CC can be mainly explained by the distance-dependent connectivity, but the common neighbor property also makes a small yet significant contribution. To illustrate this, we first investigate how the change of the decay constant of the exponential connection probability function affects the motifs; we find that when it is decreased from 12 to 4, the number of motifs with non-zero CC relative to the null model increases (r = 0.82; p < 0.001). Aside from this spatial mechanism, we find that, a stronger common neighbor effect, as captured by a common neighbor coefficient a_Γ increased from 1 to 3, also increases the number of motifs with non-zero CC value (r = 0.20; p < 0.001). These results are consistent with the observation from a recent in silico study that the overrepresented triad motifs cannot be fully captured by the distance dependence and complementary mechanisms are needed [10]; our modeling study suggests that the common neighbor property can serve as such a mechanism.

Emergence of rich-club connectivity.

We next show that our model gives rise to a rich-club connectivity phenomenon, and can explain the emergence of such a higher-order connectivity pattern. For this purpose, we calculate the normalized rich-club coefficient Φ(k), which is the ratio between the number of connections among neurons whose total degrees (i.e. the sum of in-degrees and out-degrees) are larger than k and that expected from random networks with the same in- and out-degrees for each neuron [30, 31]. In other words, this coefficient Φ(k) measures the extent to which the neurons with total degrees > k in a network are over-connected to each other. In our network, Φ(k) is consistently greater than 1.30 for neurons with k > 500 (Fig 2B blue). This result thus indicates that these neurons are over-connected to each other with at least over 30% more connections than expected from their in- and out-degrees, forming a densely connected rich-club. To further characterize these rich-club neurons, we calculate the dynamical importance for each neuron, defined as the fractional change in the largest eigenvalue of the adjacency matrix of the network upon removal of the neuron from it [32, 33], and find that rich-club neurons have higher dynamical importance (Fig 2B red). The heterogeneity of connection is also evidenced in the effectivity connectivity. To demonstrate this, we use the methods in [16, 34] to calculate the transfer entropy-based effective connectivity in our model, and find that the distribution of effective connection strengths is lognormal and there exists rich-club connectivity (see S1A and S1B Fig).

The emergence of such a rich-club connectivity structure in our model is dependent on the weakly correlated broad in- and out-degree distributions as well as the common neighbor property [3]. The former gives rise to a small group of well-connected “hub” neurons characterized by their large total degrees, and the latter makes those hub neurons more connected to each other than expected by chance. More specifically, we find that when the correlation ρ_K between in- and out-degrees or the hybrid parameter q (which controls the broadness of the in- and out-degree distributions; see Materials and methods) increases, the standard deviation of the total degree distribution also increases (r = 0.84 and r = 0.99, respectively; both p < 0.001). The common neighbor property increases the connection probability between a pair of hub neurons with large total degrees, because they are likely to share more common neighbors simply as they have more neighbors than those less well-connected neurons in the first place; indeed, we find that, as the common neighbor coefficient increases, the corresponding maximum normalized rich-club coefficient also increases (Fig 2C; r = 0.24, p < 0.001).

Rich-clubs have been found in macroscale human brain networks [31] and in the effective connectivity networks of mouse somatosensory cortex [16]. Very recently, at the synaptic level, rich-clubs have been found in a biologically constrained in silico model [10]. The rich-club behavior predicted by our circuit model is consistent with this line of research. Our model further suggests that the broad in- and out-degree distributions and the common neighbor factors are essential for the emergence of rich-clubs at the synaptic level, which could be tested in future experimental studies. Moreover, as illustrated in the following sections, the rich-club connectivity is related to heterogeneous neural dynamics as found in [12].

Transition from the asynchronous state to the propagating wave state.

By studying the spatiotemporal activity patterns emerging from our local cortical circuit model, we next illustrate that the change identified in firing rates (Fig 3A) is related to a change from an asynchronous firing state to propagating wave states. We find that when the I-E ratio ζ is large (ζ > ζ_c), the circuit does not show any structured patterns in spiking activity (Fig 3B); instead, it exhibits an asynchronous state. In this state, the globally averaged pairwise spike train correlation coefficient and the coefficient of variation of inter-spike intervals of the excitatory population are 0.002 ± 0.001 and 0.784 ± 0.002, respectively, which are similar to those found in the irregular and asynchronous firing state of randomly connected, balanced networks, e.g. [19, 36, 37].

On the other hand, when the I-E ratio ζ is small (ζ < ζ_c), coherent patterns emerge from the circuit, in the form of a localized propagating wave (Fig 3C and S1 Video; see Materials and methods for pattern detection methods), or a global plane wave if ζ is further decreased (i.e. if the circuit is further disinhibited; Fig 3D and S2 Video); these patterns propagate across the neural circuit with relatively smooth and regular moving trajectories. Such global propagating waves with regular dynamics emerging from our model due to disinhibition are largely consistent with empirical studies reporting that when cortical circuits were disinhibited by using bicuculline to block GABA receptors, similar epileptiform global plane waves were observed [38]. Such spontaneous, regular waves arising due to disinhibition are thus fundamentally different from transient propagating waves as observed in the visual cortex of behaving monkeys [39] and in the MT area of anesthetized marmosets [40, 41].

When the circuit is around the transition point between the asynchronous state and the localized propagating wave state (ζ ≈ ζ_c), spatially coherent patterns emerge intermittently from the asynchronous background activity (Fig 3E; Bayes factor > 10², see Materials and methods). Such wave patterns have complex dynamics; they propagate coherently for a while (with a mean duration of 16.5 ± 1.3 and 10% duration values larger than 45.0 ms), and then jump to a different location (Fig 3E; see S3 Video). It is interesting to note that such propagating localized patterns have been found in cortex under normal excitability conditions, and their jumping behavior as observed in our model has been explicitly documented in [13]; similar intermittent wave patterns have also been observed in the cortex of awake rabbits [42]. Note that spikes in our model are sparse (typically under 6 Hz), consistent with the spontaneous firing rate of cortical circuits in vivo [43, 44]; this thus gives rise to the propagating pattern that is spatially sparse.

We proceed to investigate the relationship between the transition state identified in our local circuit model and near-criticality dynamics. As in [45], we calculate the susceptibility χ = 〈ρ²〉 − 〈ρ〉², where 〈⋅〉 denotes averaging over time, and we use the faction of active excitatory neurons in the model for the density of activation ρ. As shown in Fig 4A, the susceptibility is the largest in the transition region from the asynchronous to the localized propagating wave state of our model. We also calculate the branching parameter σ* of spikes of excitatory neurons, which is the average ratio between the number of spikes in one bin, divided by the number of spikes in the previous bin, with different bin sizes as in [46]. It has been shown that with the bin size normalized by the averaged inter-event interval 〈IEI〉 (the inverse of population firing rate), the curves of branching parameter σ* versus the normalized bin size are very similar across species and experiments (see Fig. 7C in [46]). At the transition point of our model (as shown in Fig 4B), the curve of branching parameter σ* (subsampling 100 neurons as in [46]) is quantitatively comparable to that reported in [46]. In particular, the maximum of branching parameter σ* under subsampling in our model is 1.61 ± 0.06, which is close to that of 1.4 in [46]. These results thus indicate that the transition state in our model is near criticality.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Near critical dynamics around the transition of the circuit activity states.

(A) Susceptibility of the excitatory population at different I-E ratios ζ. 1 ms time bin is used for the calculation; the results are not sensitive to the choice of bin size. The error bars show one SEM of 10 trials. (B) Branching parameter σ* under subsampling as a function of the normalized bin size at ξ/ξ_c = 1. The shaded area shows one SEM of 10 trials.

https://doi.org/10.1371/journal.pcbi.1006902.g004

Theoretical analysis.

To understand the mechanisms underlying this transition from the asynchronous state to the localized propagating wave state, we employ an analytical framework known as the diffusion approximation [47, 48], in particular, its recent extension for analyzing spatially extended, spiking neural circuits [49] (see S1 Appendix for details). Unlike more abstract models commonly used for studying emergent spatial patterns in neural media, this analytical approach preserves most of the biophysical details of the spiking neurons (Eq 4). We find that when the mean E-to-E synaptic strength used in the analysis is scaled by an empirical factor h = 0.4775, the critical I-E ratio ζ_c from the analysis and the simulation coincide; the role of h is to compensate for the simplifications imposed by the diffusion approximation (see S1 Appendix for details). As shown in Fig 5A, for ζ > ζ_c, the mean firing rate of the spatially uniform activity predicted by the analysis, which corresponds to the asynchronous state, matches closely to our full spiking circuit model; however, for ζ < ζ_c, the circuit’s mean firing rate becomes larger than that predicted by the analysis, as a sharp transition to spatially coherent pattern occurs in the full circuit model. We analyze the stability of the spatially uniform activity under spatially periodic perturbations, and find that it is stable for ζ > ζ_c but unstable for ζ < ζ_c. In the latter case, we find that the formation of spatially coherent pattern is due to a Turing bifurcation (Fig 5B). Our analysis thus provides an explanation for the formation of the spatially coherent pattern in our cortical circuit model.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Stability analysis of the homogeneous, asynchronous state of the spatially-extended spiking circuit.

(A) Mean firing rate of the excitatory population in the simulation of our local cortical circuits (dots) and the rate of spatially uniform activity from the analysis (solid line) with the optimal correction factor. The dashed line shows the critical point ζ = ζ_c. (B) Eigenvalues with the largest real parts of the network’s dynamics in response to spatially periodic perturbations (spatial Fourier mode k = (1, 1)), with different I-E ratios ζ and the optimal correction factor. The eigenvalues in this parameter regime are real, with positive values emerging as the I-E ratio ζ decreases below the critical point ζ_c (dashed line), indicating a Turing bifurcation. The dominant spatial eigenmode has wave vector k = (1, 1), indicating that a single pattern is formed, consistent with that observed in the full spiking circuit.

https://doi.org/10.1371/journal.pcbi.1006902.g005

Common neighbor property of neural connectivity can induce the change of dynamic states.

We now demonstrate that the common neighbor factor provides a novel structure-based mechanism for shifting cortical activity states. We vary the common neighbor coefficient a_Γ; a larger a_Γ means that connections will more likely be formed between neuron pairs sharing more common pre-synaptic neighbors (see Materials and methods). Note that the value of a_Γ does not affect the number of connections or the total incoming connection strengths of individual neurons. Instead, a larger a_Γ results in a more clustered network (i.e. with more triad motifs that have closed loops), as indicated by an increasing mean network CC (Fig 1D). We find that such a change in the connectivity induces a similar transition from the asynchronous state to the localized propagating wave state. Networks with smaller a_Γ and thus smaller CC values fire asynchronously and irregularly (Fig 6A); on the other hand, networks with larger a_Γ and CC values exhibit localized propagating patterns with regular moving trajectories (Fig 6C). Propagating patterns with complex dynamics, including variable propagation direction, speed, and jumping behavior, thus emerge from the network with intermediate CC values (Fig 6B); this scenario is largely consistent with the recent in vitro finding that moderately clustered networks appear optimal for initiation and propagation of diverse, complex propagating wave patterns [50].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Common neighbor coefficient a_Γ can induce the transition from the asynchronous state to the localized propagating wave state.

(A-C) Snapshots of spiking patterns in the local cortical circuit model with different common neighbor coefficient a_Γ values. These patterns are visualized in the same way as in Fig 3B–3E. (D) Phase diagram constructed based on the pattern detection rate (the percentage of time frames where a localized wave is detected) and the average firing rate of the excitatory population. The States 1, 2, 3 and 4 correspond to the asynchronous state, the transition state, the localized propagating wave state and the global plane wave state, respectively. The circle denotes the working point of our local cortical circuit model. (E) Phase diagram as in (D) but with the spike frequency adaptation removed, where the different states are shown by red boundary lines. The state I, II, and III correspond to the asynchronous state, the state with localized activity patterns that barely propagate (see S4 Video), and the global plane wave state, respectively. The excitatory neurons without firing a spike during the 10 second simulation time are detected as quiescent neurons.

https://doi.org/10.1371/journal.pcbi.1006902.g006

These results reveal that both the neurophysiological mechanism (i.e. the I-E ratio ζ) and the structure-based mechanism (i.e. the common neighbor coefficient a_Γ) can shift circuit activity state from the asynchronous to the propagating wave state. We thus use these two parameters (ζ and a_Γ) to systematically explore the circuit dynamics based on the distinctive dynamic states we have found above. As shown in the phase diagram (Fig 6D), the transition from the asynchronous to localized propagating wave state can be induced across a range of the parameter values. We define a working point of the cortical circuit model at one such critical transition of the states with ζ = ζ_c and a_Γ = 2 (Fig 6D), where the average firing rate (3.0 ± 0.2) of the excitatory neurons are close to that found during spontaneous activity in vivo [43, 44].

Effect of spike-frequency adaptation.

We now demonstrate that the spike-frequency adaptation (Eq 5) affects the spatiotemporal dynamics of our local cortical circuit model. For this purpose, we investigate neural circuit dynamics after removing the spike-frequency adaptation from our model (by setting Δg_k = 0 in Eq 5). As shown in Fig 6E, spatiotemporal patterns can still be formed without adaptation, including localized activity patterns (State II) and global plane waves (State III). However, the localized activity patterns in State II have different dynamics, which are pinned to certain locations without clear long-range propagations as found in the original model with adaptation. This behavior of the localized patterns results in the scenario that within this state, only a small portion of excitatory neurons fire and the rest is quiescent. As in [9], we calculate the fraction of excitatory neurons quiescent in our spatially extended, heterogeneous network and find that this fraction is mostly greater than 50% for this state. In contrast, for the original network with the adaptation, the fraction of neurons quiescent is low (< 0.1%) for all the states. This result is consistent the observation that the inclusion of adaptation can reduce the fraction of neurons quiescent in randomly connected neural circuits with heterogeneous degree distributions [9]. In comparison, the removal of adaptation has less effect on the asynchronous state (State I in Fig 6E); the part of the asynchronous state in the original phase space with lower firing rates (Fig 6D; bottom right), after the removal of adaptation, remains as the asynchronous state with dynamics similar to those found in [37]. This is expected because in this region, the inter-spike intervals of individual neurons are rarely shorter than or comparable to the decay constant of potassium current (τ_K = 80 ms; see Eq 5), rendering the effect of spike-frequency adaptation negligible. Note that both an increase in the I-E ratio ζ in the circuit without spike-frequency adaptation and the introduction of adaptation can induce the transition from the state with localized activity patterns that barely propagate to the asynchronous state; this suggests that both factors tend to reduce excitation in the circuit. However, only the introduction of adaptation can induce the transitions from local to propagating patterns. This happens because the adaptation destabilizes the local pattern and fluctuations of the network would cause it to propagate away from its original location, forming a propagating wave.

Emergent, tightly balanced excitation and inhibition.

It has been found in vivo that cortical neurons receive strongly correlated excitatory and inhibitory inputs, with the former closely tracked by the latter as characterized by a strong cross-correlation around zero time lag [23, 51]; such tight proportionality between E and I inputs is regarded as a tight balance [52]. We now show that tightly balanced synaptic inputs with the dynamical properties as found in [23] emerge from the working point of our model as defined above; the balanced condition in our current model is thus referred to as such a tightly balanced one. Fig 7A (middle panel) shows that the magnitudes of E and I currents into a neuron closely track each other in time, with strong cross-correlations around zero time lag (Fig 7B red), which decay to zero within around 50 ms; note that this decay time is largely consistent with that measured in in vivo recordings (see Fig. 1 of [23]). It is also interesting to note that only near the critical transition state, as shown in Fig 7A (middle panel), synaptic inputs appear to exhibit large, occasional excursions in amplitudes; such bumpy properties of synaptic inputs have been explicitly documented in [23]. When the system is in the asynchronous state, the inputs do not have such bumpy properties (Fig 7A bottom panel) but the E and I currents are still highly correlated at zero time lag. However, the cross-correlation decays to zero much faster (within around 15 ms; Fig 7B blue) than that measured in [23]. In the propagating wave state, the strong cross-correlation is also present (Fig 7B cyan) but the synaptic inputs become semi-periodic and regular (Fig 7A top panel), apparently inconsistent with the large fluctuations observed in synaptic inputs in vivo.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Properties of synaptic inputs to individual neurons in the local cortical circuit model.

(A) Time series of total excitatory and inhibitory inputs received by individual excitatory neurons over 1 s in the different dynamic states (with ζ/ζ_c = 0.8, 1.0, 1.3 corresponding to the propagating wave state, the transition state and the asynchronous state, respectively; see Fig 3). (B) Average cross-correlation (xcorr) between the total E and I currents into excitatory neurons in the different dynamic states. The shaded area shows one SEM of 10 trials. (C) The convex hulls of the mean and standard deviation data points (pooled from 10 trials) of the total input current into individual neurons in the different dynamic states, denoted by the same color-coding as in (B). The region to the left of vertical dashed line is where the mean input current is less than the threshold current. The region above the inclined dotted line is where the mean plus one standard deviation of the input current is above the threshold current. (D) The relative errors of estimating the total standard deviation σ_I (and mean μ_I) of the recurrent currents with the sum of individual components as given by the variance decomposition (Variance decomposition analysis), averaged over the excitatory neurons and over a period of 50 ms within each trial. (E) The contributions to the variance in recurrent currents from various sources are shown across circuits with different I-E ratio ζ. The numbers above the bars are the corresponding total variances averaged over excitatory neurons (nA²).

https://doi.org/10.1371/journal.pcbi.1006902.g007

Another key property of the synaptic inputs to cortical neurons is that they keep neurons in a high-conductance state [24, 53]. We now demonstrate that in our conductance-based model, the statistical properties of such a high-conductance state can be captured. At the working point, the mean excitatory synaptic conductance into individual excitatory neurons is of the order of the leaky conductance (with the ratio between them being 1.27 ± 0.15; 400 sample neurons pooled from 10 trials), as found in pyramidal neurons in areas 5, 7, and 21 of the parietal cortex of awake cats [24]. Furthermore, the mean and fluctuation (i.e., SD) ratios between inhibitory and excitatory conductances (2.81 ± 0.21 and 2.95 ± 0.59, respectively) are also consistent with those measured in [24]. As a result of such synaptic inputs, the membrane potential is depolarized (with a mean of −61.7 ± 0.6 mV) and exhibits large fluctuations (with a SD of 4.1 ± 0.3 mV). However, our model does not capture the observed decrease in inhibitory conductance within 20 ms before a spike [24]. As the I-E ratio becomes larger (ζ > ζ_c), and mean excitatory conductance is still comparable to the leaky conductance (with a ratio of 0.79 ± 0.03 at ζ/ζ_c = 1.3, in the asynchronous state); however, when the circuit becomes disinhibited with a smaller I-E ratio (ζ < ζ_c), the mean excitatory conductance becomes pathologically high (6.67 ± 0.68 times the leaky conductance at ζ/ζ_c = 0.8, in the propagating wave state).

Structural heterogeneities contribute to neural fluctuations.

Similar to the conductances, the excitatory and inhibitory input currents into individual neurons also show large fluctuations (Fig 7A). This dynamic regime near the working point of our model is fluctuation-driven [54]; that is, the mean input currents into individual excitatory neurons are subthreshold as a result of the tight balance described above and, without the time varying fluctuations around their mean input, the neurons would not fire (Fig 7C red). In this regime, the neural dynamics are largely determined by the variances in the recurrent currents and therefore highly variable, as characterized by a mean coefficient of variation of inter-spike intervals of the excitatory neurons at 0.93 ± 0.03, consistent with that measured from the spontaneous activity in vivo. In the asynchronous state (ζ > ζ_c), the mean input currents remain subthreshold (Fig 7C blue); on the other hand, in the disinhibited propagating wave state (ζ < ζ_c), the mean input currents become superthreshold (Fig 7C cyan).

By incorporating a range of realistic features of cortical connectivity into the spiking circuit model, such as the heterogeneities in both connection topology and strength, our model allows us to explore their effects on the variable neural dynamics. In particular, we investigate how the structural heterogeneities contribute to the variances in the recurrent excitatory currents in the fluctuation-driven circuits with ζ ≥ ζ_c. For this purpose, we develop an analysis method called variance decomposition analysis (see Materials and methods). As described in Eq 9 syn, there are four components determining the recurrent currents, namely, the connection topology, the connection strength, the pre-synaptic dynamics and the postsynaptic dynamics, denoted by the shorthand letters a, J, s and V, respectively. Our analysis shows that there are 15 different sources of variance of the recurrent currents (Variance decomposition analysis), arising from each of the four components (a, J, s and V) contributing via either its own (co-)variance or mean value (see Materials and methods). For example, one such source arises from the coupling among the variances of s & a and the mean values of V & J; we denote this source by , only including in the subscript the letters representing those components contributing via their (co-)variances. In our variance decomposition analysis we assume that the four components are mutually independent and ignore the moments beyond second-order. To validate this analysis, we study how well the sum of individual variance sources (RHS of Variance decomposition analysis) matches the total variance (LHS of Variance decomposition analysis) using numerical data with different I-E ratio ζ. As shown in Fig 7D, the two match with reasonable accuracy.

The contributions of these 15 sources across fluctuation-driven systems with different I-E ratios are shown in Fig 7E; we find that, to capture most of the variance in recurrent currents across different I-E ratios, at least four sources need to be considered, including , , and . Note that the variances in connection topology (a) and strength (J) contribute to 3 out of 4 of the major sources; as shown in Fig 7E, without these two structural factors, at ζ = ζ_c, only 78.3 ± 9.9% of the total variance can be captured by alone. Therefore, the variances of connection topology and strength make a significant contribution to the variable recurrent currents at the working point of our model. When the I-E ratio increases to ζ/ζ_c = 1.4, the amount of variance captured by alone further decreases down to around 50%. This result thus reveals the importance of the structural heterogeneities in generating the variable, fluctuation-driven neural dynamics, going beyond previous studies about variable neural dynamics, in which such variances in connection topology and strength were not captured, e.g. [19, 36, 47].

Emergent heterogeneous population coupling.

Recently, it has been found that in the deep layers of sensory cortex of awake mammals, individual neurons correlate (couple) to the population activity differently, ranging from strongly correlated “choristers” to weakly correlated “soloists” [12]; this experimental study has revealed a fundamental relationship between individual neurons and a larger population. However, the dynamic nature of this diverse coupling behavior and how it is related to structural connectivity remains unclear. We now demonstrate that such diverse coupling can emerge in our model near the critical transition from the asynchronous to the localized propagating wave state. In addition, we illustrate that these heterogeneous neural dynamics are related to the rich-club connectivity predicted by our model.

To quantify the population coupling in our model, we calculate the spike-trigged population rate (stPR); as in [12], this is done by first summing the spike train of mutiple neurons to get the population rate, and then correlating the spike train of each neuron with this population rate. As shown in Fig 8A, the neurons in our model exhibit diverse population coupling behavior; that is, there are neurons with strongly positive, weakly positive, neutral or even negative correlations with the population activity. To test the significance of this diverse coupling, we use the method in [12] to shuffle the spike data in a manner that preserves both the population rate and the average firing rate of each neuron but randomizes the timing of individual spikes. We find that, after the shuffling, the diverse population coupling behavior almost completely disappears (Fig 8B). For a quantitative comparison, we calculate the stPR coefficients as in [12]. In our model, the stPR coefficients before shuffling ranges from −0.96 ± 0.26 to 2.85 ± 0.23, with a mean of 0.81 ± 0.09 and a standard deviation of 0.76 ± 0.06; these results are comparable to those measured in the deep layers of S1 of awake mice, i.e. with a range from -0.5 to 2.5 and a mean around 0.5 [12]. After shuffling, however, the standard deviation of stPR coefficients becomes much smaller (0.14 ± 0.03; Fig 8B inset; p < 0.001, two-sample F-test), confirming that the diverse population coupling behavior emerging from our model is an inherent dynamical property of our local cortical circuit model.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 8. Emergent, heterogeneous population coupling in the dynamical working regime of the local cortical circuit model.

(A) Spike-triggered population rate (stPR) for five representative neurons in a single trial. (B) stPR for the same five neurons but with spikes shuffled. Inset: distributions of stPR coefficients pooled from 10 trials before (black) and after (red) spike shuffling. The coefficients are defined as the stPR values at zero time lag normalized by the median value of the shuffled data. For each trial, 10 samples of 66 excitatory neurons are randomly chosen to match the sample size used in experimental studies [12]. Neurons are sampled from a circular area with a radius of 25 grid units (about 125 ¼m) for approximating the spike detection range of a single shank of silicon electrode array used in [12]. The results are not sensitive to the size of the sampling area or the number of sampled neurons. (C-D) stPRs for ζ = 1.15 and ζ = 0.7, respectively. Insets: shuffled stPRs. For ζ = 1.15, the model is shifted away from the working regime into the regular propagating wave state; while for ζ = 0.7, it is shifted into the irregular, asynchronous state.

https://doi.org/10.1371/journal.pcbi.1006902.g008

We find that our model does not show such diverse population coupling when it is shifted away from the working point (i.e., the critical transition point between the asynchronous and the propagating wave states), although the structural topology remains the same. As shown in Fig 8C, in the asynchronous state, the stPRs are noisy and show little difference from the shuffled results. As the network becomes more disinhibited with a smaller ζ value, the stPRs become periodic due to the regular movement of patterns (Fig 8D), not consistent with those reported in [12].

To unravel whether and how the heterogeneous connectivity is related to this diverse coupling, we compare the rich-club connectivity of the neurons with their population coupling behavior. As mentioned above, the normalized rich-club coefficient Φ(k) characterizes the degree of “over-connection” among the neurons with a total degree higher than k; Fig 9A shows the connections among such a group of high degree neurons. In order to compare Φ(k) against the population coupling behavior, we need to quantify the corresponding degree of “over-coupling” to the population rate of the corresponding neurons (with a total degree higher than k). This can be achieved by calculating the average z-score of stPR coefficients for those neurons. As shown in Fig 9B, this k-dependent stPR coefficient z-score shows an overall increasing trend when the total degree k increases; this trend is closely followed by the normalized rich-club coefficient Φ(k). For a more direct comparison, we plot the coefficient z-scores against the corresponding normalized rich-club coefficients with the same k (Fig 9B inset), and confirm that the z-scores can be largely explained by Φ(k) (r = 0.84; p < 0.001).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 9. Rich-club connectivity is related to the heterogeneous population coupling.

(A) An example of the rich-club connectivity among high total degree (k > 600) excitatory neurons (red dots). Only a 10% random subset of the synaptic connections are plotted (lines) for visualization purpose, with a probability of 1/Φ(k = 600) colored in blue and otherwise in yellow to illustrate the proportion between the number of connections expected from their total degrees alone (blue) and the number of extra connections due to the rich-club connectivity (green). The circle shows the range within which both the rich-club and stPR coefficients are calculated as used in Figs 2B and 8. (B) The rich-club neurons have high stPR coefficients. Left y-axis (blue): the blue dashed line shows the stPR coefficient z scores for the neurons with a total degree > k averaged over 10 trials and the shaded area shows one SEM. Within each trial, 10 random 66-neuron samples (as in Fig 8) are collected to calculate the k-dependent stPR coefficient z scores. Right y-axis (red): the red dashed line shows the normalized rich-club coefficient curve from Fig 2B and the shaded area shows one SEM. Inset: normalized rich-club coefficient Φ(k) (y-axis) versus k-dependent stPR coefficient z score (x-axis) from 10 realizations of the network (shown with 10 different colors). (C) Distributions of stPR coefficients pooled from 10 trials before (black) and after (red) spike shuffling from the local cortical circuits with different common neighbour factors a_Γ. For comparison, the total count (area under curve) and bin size are the same for each distribution. (D) stPR coefficients before shuffling, normalized by the SD of stPR coefficients after shuffling (i.e., coefficients of control case) and plotted in the same k-dependent manner as in (B). Each curve shows the average of 10 trials.

https://doi.org/10.1371/journal.pcbi.1006902.g009

To further illustrate the relationship between the rich-club connectivity and the diverse population coupling, we vary the common neighbour factor a_Γ that can modulate the rich-club coefficients (Fig 2C) while preserving the in- and out-degrees of individual neurons (see Materials and methods). We find that when the common neighbour factor a_Γ decreases from 3 to 1, the difference between the original distribution of stPR coefficients and the distribution of the control case (i.e., the case with random shuffling of spikes) gradually diminishes (Fig 9C), indicating that the diverse population coupling effect becomes less significant. This trend can be quantified by normalizing the original stPR coefficients with the standard deviation (SD) of those coefficients of the control case. These normalized coefficients are then compared across circuits with different values of a_Γ; as shown in Fig 9D, the trial-average of such coefficients decreases as the common neighbour factor a_Γ decreases. As decreasing a_Γ results in a decrease in the peak value of normalized rich-club coefficient Φ(k) (see Fig 2C), these results indicate that the rich-club connectivity in our circuit model is closely related to the diverse coupling behavior of neural dynamics.

Dynamic propagating patterns and precise spiking structures.

We next demonstrate that in our local cortical circuit model, propagating wave patterns emerging near the transition point (Fig 6) can quantitatively account for the dynamic properties of propagating spiking waves of UP states found in spontaneous activities of cortical neurons in layer 5 of awake rats [14]. In addition, our model can capture precise spiking structures, which have been found to coexist with propagating waves [14]. We further show that the spatial extension property of our model is essential for the emergence of these patterns.

Around the working point (i.e. the transition point from the asynchronous to the localized propagating wave state), our model exhibits a propagating wave pattern with irregular trajectory and variable speed as illustrated above (Fig 3E). To compare these wave patterns with those found in [14], we visualize the firing activity using raster plots where neurons are sorted based on one spatial coordinate; Fig 10A shows a wave propagating with a component of its direction along this spatial coordinate. To characterize the pattern propagation in our model, we calculate the distance traveled every millisecond to obtain the propagation speed distribution. We find the propagation speed is highly variable, as evidenced by its log-normal distribution (Fig 10B). This propagation speed distribution has a mean of 1.72 ± 0.06 grid-units ms⁻¹ (12.7 μm ms⁻¹; see Materials and methods for the conversion) and 10% of values larger than 3.55 ± 0.11 grid-units ms⁻¹ (26.3 μm ms⁻¹), which matches the highly variable propagation speed measured in the cortex of awake rats with values spanning from 8 to 40 μm ms⁻¹ [14]. Such variability in speed has also been found in cortical slices, where the propagation of localized excitation patterns often abruptly change speed; under normal excitability conditions, the typical speed values are between 6-10 μm ms⁻¹ [13]. In our model, the variability of propagating wave patterns is lost when the system is disinhibited (with small I-E ratio ζ), as patterns propagate in a more regular way.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 10. Dynamic properties of propagating wave patterns emerging in the dynamical working regime of the local cortical circuit model.

(A) Raster plots of 50 neurons randomly sampled from a 10-by-50 grid-unit region, sorted according to their y-coordinates, showing two travelling patterns moving in opposite directions at different time points. The solid blue line is the linear fit to the first spikes at each y-position (y = 0 denoting the y-center of the sampling region). (B) The distribution of the pattern propagation speed (dots) with the log-normal fit (blue solid line). The error bars and shaded area show one SEM.

https://doi.org/10.1371/journal.pcbi.1006902.g010

Around the working point, our model can also capture the precise spiking structures during spontaneous activity in layer 5 of rodent somatosensory cortex [14]. Precisely repeating patterns embeded in variable neural activity have also been found in spontaneous activity with critical dynamics [55]. In our model, precise spiking structures can be detected by considering spike triplets. As illustrated in Fig 11A, a spike triplet is composed of three spikes, each from a different neuron, and can be characterized by the two inter-spike intervals. The joint distribution of the two inter-spike intervals should have a flat structure without any preferred mode if the three neurons fire randomly as independent stationary Poisson processes. Therefore, a statistically unlikely large peak in the joint distribution, as found in our model (Fig 11B), means a temporally precise, sequential firing pattern appears repeatedly in this neuron trio. Following the method in [14], the triplets with inter-spike intervals 10 ms around the mode are detected as precisely repeating triplets.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 11. Emergent precise spiking structures.

(A) Schematic of a spike triplet that is described by two inter-spike intervals. (B) The probability distribution of two inter-spike intervals from a representative neuron trio. The black square denotes the precisely repeating triplets occurring around ± 10 ms of the mode. (C) Occurrence of precisely repeating triplets peaks shortly after the start of activated states. The time from activated state (AS) onset is calculated based on the first spike of the triplet. The black, red, and green lines are the trial-averaged normalized probability density functions (PDFs) for the original data, shuffled data [14], and data dithered with different window sizes (method L in [56]), respectively. For each trial, 50 samples of 50 excitatory neurons are randomly chosen from a 10-by-50 grid-unit region (about 50-by-250 μm) to approximate the spike detection region of two neighboring shanks of silicon electrode array used in the experimental study by [14]. The results are not sensitive to the width of the rectangular sampling area. Within each sample, 50 neuron trios are randomly chosen to get the original and shuffled/dithered spike triplet counts.

https://doi.org/10.1371/journal.pcbi.1006902.g011

To further test the statistical significance of the detected precisely repeating triplets in our model, we first shuffle the individual spike timing using the method in [14] and compare the number of detected triplets before and after the shuffling, with respect to the time from activated state onset, which serves as an objective temporal reference for comparison. As shown in Fig 11C, the number of the detected triplets (or equivalently the triplet probability) is consistently lower after shuffling within 100 ms after the activated state onset, confirming that a significant portion of detected triplets in the original data are intrinsically structured. In particular, the triplet probability before shuffling peaks shortly after the activated state onset (around 48 ms; peak value 0.78 ± 0.50 times higher than that of the shuffled data; p < 0.001 for two sample t-test), which is consistent with the in vivo findings [14]. Additionally, we use the dithering method (method L in [56]) to test the statistical significance of this result. As shown in Fig 11C, the dithering method gives a characteristic decay of the triplet probability when the dithering window size increases from 10 ms to 80 ms (compared with the triplet detection window size 20 ms), and the peak triplet probability before dithering is significantly higher than that of the dithered ones (p < 0.006 for two sample t-test), confirming the statistical significance of the detected triplets [57].

The detected activated states largely correspond to a spiking wave pattern passing through the sampling region, as the chance of finding the center of a wave pattern within the sampling region during activated states is nearly three times (2.84 ± 0.62) as high as the average. Therefore, the results in Fig 11C indicate that the precisely repeating triplets are initiated by the propagating waves patterns, as found in [14]. In our model, when the I-E ratio ζ is shifted further into either the asynchronous state or the regular propagating wave state, no such statistically significant precise spiking structure can be detected. In [58], however, spike sequences can reliably arise in an irregular, asynchronous state of randomly connected neural networks with lognormal synaptic weight distributions, and these sequences are organized along the unidirectional synaptic pathways embedded by rare strong synapses.

In our local cortical circuit model, due to its spatial extension property formed by the distance-dependent connectivity, spikes of one group of neurons would result in those neurons that are physically proximal to them to fire, resulting in propagating patterns. To further dissect this role of spatial extension in our model, we set all of the decay constants of the exponential connection probability functions in our model , so that the connectivity no longer depends on the spatial locations of the neurons. This effectively removes the spatial extension property, while keeping all the other connectivity properties of our model fixed. We find that propagating patterns cannot be formed and precisely repeating triplets cannot be detected. These results thus demonstrate the importance of such a spatial extension arising from the distance-dependent connectivity in shaping complex cortical dynamics, indicating that structured spatiotemporal patterns can emerge from the highly variable firing activity of individual neurons in such spatially extended networks.

Heterogeneous degree distributions.

We now specify the heterogeneous degree distributions used in Step 2 of our generative connectivity model. Direct evidence of heterogeneous degree distributions of neural connectivity at the synaptic level is lacking, mainly due to small sample sizes as limited by currently available tracing techniques [25]. However, heavily-tailed degree distributions have been found in the effective connectivity of rodent somatosensory cortex using transfer entropy-based methods [7], which produce qualitatively similar results to those found by the patch clamp studies of [3]. Such degree distributions have also been observed in a biologically constrained in silico microcircuit of the same area [10]. Anatomically based estimates of the connectivity of layer 4 rat barrel cortex [9] have also indicated that the degree distributions of cortical neurons are significantly broader than Poisson distributions in the conventional E-R random networks (where the connection probabilities between any pair of neurons are identical) for a given mean value. In light of these studies, in our model, we assume that both the in- and out-degrees are heterogeneous, with some neurons having considerably more incoming and/or outgoing connections than expected in the E-R networks. More specifically, in our model, the in-degree K_in and out-degree K_out for any excitatory neuron are sampled from a hybrid distribution composed of a linear combination of a Poisson distribution and a log-normal distribution, which is a mathematically simple and well-known heavy-tailed distribution, (3) where K_in and K_out are rounded to the nearest integers, q is the hybrid parameter, and are drawn from Poisson distributions with mean ζ = p^EEN^E, p^EE = 0.16 is the overall E-to-E connection probability; and are drawn from log-normal distributions with the same mean ζ and a standard deviation that is 20% of the mean. Furthermore, to allow for a control of the correlation between in- and out-degree as a system parameter, the two Poisson random variables and the two log-normal random variables are generated with identical correlation coefficients of ρ_K; in our model, we assume there is a weak positive correlation ρ_K = 0.13; note that positively correlated in- and out-degrees have been found in microscopic connectivity data [25]. In our model, the hybrid parameter q allows the interpolation between Poisson distributions as in E-R networks (q = 0) and log-normal distributions (q = 1); with the default hybrid parameter q = 0.4, the hybrid in- and out-degree distributions have a mean and standard deviation of 635 and 83, respectively. The total sums of in- and out-degrees in a network must be exactly equal, which is not guaranteed by sampling them from distributions with the same means; to satisfy this constraint, a random subset of sampled in- or out-degrees are adjusted by one so that any mismatch in the total sums is compensated for.

Equalized I-E ratios.

In our model, we consider an essential neurophysiological feature of local cortical circuits, that is, the excitatory post-synaptic currents are proportional to the inhibitory ones, with a homogeneous ratio across all excitatory neurons, as found in mouse primary visual cortex layer 2/3 [18]. To model this in our heterogeneous circuit, we consider the I-E ratio , where denotes the number of connections (in-degree) received by excitatory neuron i from the inhibitory population and the connection strength values are determined by the reverse pooling method mentioned above. To equalize the I-E ratio ζ_i across the neurons to a desired network-wide ratio, that is, 〈ζ_i〉 = ζ, the values for neuron i are sampled from a Gaussian distribution with a mean equal to and a standard deviation that is 25% of the mean. The I-E ratio ζ will be varied as a system parameter to explore the spatiotemporal dynamics of our cortical circuit model.

We now summarize the parameters used for generating the entire network. For most of the results, we use N^E = 63 × 63 = 3969 excitatory neurons and N^I = 1000 inhibitory neurons. The ratio between N^E and N^I is around 4, consistent with published measurements from rat visual cortex [73]. We consider the distance between two neighboring excitatory neurons (one grid unit) in our model to be 7.4 μm based on the measurements from rat primary visual cortex layer 5 [74]. Therefore, our network model contains approximately the same number of pyramidal neurons as a 0.22 mm² square patch of cortical layer 5. Our results are not sensitive to the size of the network; we have obtained similar results for larger networks with sizes such as N^E = 103 × 103 = 10609 and N^I = 2673. For the exponential distance-dependent connection probability functions used in Step 1 of our generative connectivity model, the decay constants for the excitatory connections are grid units and grid units (corresponding to roughly 74 μm), which are consistent with experimentally measured ranges [3]. For the inhibitory connection ranges, we use grid units. For simplicity, the in-degree distributions for the connections from and to the inhibitory population are Poisson with overall connection probabilities p^EI = p^IE = 0.2 and p^II = 0.4, taking into account the fact that the inhibitory connections are locally denser than excitatory ones [75, 76]. The rest of the connection strengths are assumed to be homogeneous and we use and . Our connectivity model can be regarded as a generic local circuit model of mammalian cortex. However, considering the variations across species and cortical areas in general, e.g. [77], our model is most applicable to the layer 5 of rodent primary sensory cortices, whose connectivity properties are used to constrain the majority of the connectivity features in our model.

Variance decomposition analysis.

To study how different heterogeneities in the neural circuit model each impact neural dynamics, we develop the following variance decomposition analysis. As described in Eqs (6 and 7), the recurrent synaptic current in the conductance-based model is given by (9) We obtain the deviation parts of the connection topology (denoted by ) and strength (denoted by ) by removing the mean values (over different realizations), that is, and . We denote the standard deviations as and , respectively. Similarly, assuming the dynamics are stationary in terms of firing rate of each neuron, we have the time varying parts of the gating variable and the membrane potential as and , where the averages are taken over an ensemble of networks with different realizations of the random initial conditions. Following the same notation, the corresponding standard deviations are and .

Substituting the above expressions into Eq 9, expanding the brackets and after some derivations utilizing or assuming the independencies among the deviation and time varying parts (see S3 Appendix), we obtain the mean and time varying part of the recurrent current, , where the mean is given by and the time varying part can be expressed as a sum of the products of 15 random variables and their associated standard deviations, (10) where are independent normal random variables with zero mean and unit variance and the associated standard deviations can be systematically obtained and expressed in terms of the mean or (co-)variance of the four components in Eq 9, namely, and . The subscript X contains the shorthand letters of the components contributing to the term via their (co-)variances instead of mean values; for instance, for X = as, (see S3 Appendix for detailed derivations). And there are in total 15 different terms for the recurrent current from population β to α, that is, X ∈ {a, J, s, V, aJ, as, aV, Js, JV, sV, JsV, asV, aJV, aJs, aJsV}. These results enable us to obtain the variance σ²(⋅) of the synaptic currents, (11) which allows us to decompose the contributions to the variance of recurrent currents from different sources, including connectivities.

Localized propagating wave detection.

The spontaneous activity of our local cortical circuit model exhibits rich spatiotemporal dynamics, including propagating wave patterns. To compare the properties of these wave patterns in our model with those reported in [14], we use the following detection method. We assume at any time the excitatory population has a firing rate profile of a 2-D circular Gaussian function, with the firing rate of excitatory neuron i described as (12) where r_i is the coordinate vector of neuron i and is the parameter vector that contains the center position, the variance, and the height of the Gaussian function, respectively. We further assume that the firing process of each neuron around time t can be approximated by a Poisson process with rate . With this assumption, the log-likelihood ln L(θ_t) of the parameters is given by (13) where is the spike count of excitatory neuron i within a small window Δt = 5 ms centered around time t. To find the maximal likelihood fit to the spiking data on the grid with periodic boundaries, we transform the grid coordinates into points on unit circles and calculate the circular means. To measure the goodness of fit, we compare the Gaussian model against a null model that assumes a uniform firing rate profile by calculating the Bayes factor B₁₂ (see S4 Appendix for details). A Bayes factor log₁₀(B₁₂) > 2 means the Gaussian model is significantly more strongly supported by the data than the null model is [86]. In our analysis, we regard any fitting results with log₁₀(B₁₂) < 2 as irregular firing states between coherent spiking wave patterns.