Intrinsic neural diversity quenches the dynamic volatility of neural networks

We further characterized the influence of heterogeneity on the dynamics and volatility of our network in Fig. 2. In the homogeneous case (i.e., σ_H² = 0), neuronal activity was indeed characterized by alternating epochs of stability and instability, as portrayed in Fig. 2A. Transient stable dynamics indexes states in which neural activity is stable and relaxes back to equilibrium after perturbations. In contrast, other periods were characterized by unstable neural activity in which the activity of individual neurons diverges away from equilibrium. Such periods of instability result from critical transitions (15) in which neural activity departs from stability, and may diverge, become synchronous and/or chaotic.

To illustrate the corresponding activity of the neurons, we generated exemplar Poisson spike trains whose rates match those predicted by our network model (Materials and Methods). As illustrated in Fig. 2A, modulatory input triggered abrupt and frequent transitions between regimes of quiescence (i.e., low firing rates) and high-frequency firing in which burst-like activity could also be observed. Such periods of intense spiking repetitively followed the onset of instability, resulting in network-wide correlated fluctuations. The presence of such covariability is also apparent through the high variance of the mean network firing rate. These observations suggest that slowly fluctuating modulatory input may expose the volatility of homogeneous networks by revealing sudden stability transitions and dynamical regimes that are qualitatively distinct, exemplifying nonresilient behavior.

In contrast, when excitability heterogeneity is introduced (i.e., σH2 > 0) in Fig. 2B, significantly different dynamics are observed despite identical connectivity and modulatory drive. Indeed, the long-term dynamics of the network are resilient: Robust, invariant stability replaced the intermittent behavior seen in the homogeneous case. No transitions between stability and/or instability occurred. As a direct consequence of heterogeneity, redundancy in the neurons’ equilibria is broken: Neuron fixed points were now distributed with a mean μuo = B + So and variance σuo2 (Materials and Methods). As can be seen in the representative spiking activity of neurons, such differences in excitability yield increasingly moderate responses to modulation. Neurons with high excitability compensate for those with low excitability, and vice versa, thereby restraining network fluctuations by distributing them across the population. This translates into more gradual changes in the network mean firing rate, compared to the homogeneous case (cf. Fig. 2 A and B).

We quantified both numerically and mathematically such changes in mean network firing rate, to quantify network-wide excitability (Materials and Methods). Fig. 2C shows that excitability heterogeneity promotes increasingly linear, moderate changes in mean firing rate resulting from modulatory input. In contrast, homogeneous networks exhibit steep, nonlinear responses. Such results confirm that differences in excitability restrain firing rate fluctuations across the population by limiting sudden increases in activity. Counterintuitively, a consequence of these results is that heterogeneity also enhances the net excitability of the network: Smaller modulatory amplitudes are required to trigger a (albeit weaker) response compared to homogeneous networks. Conversely, decreased heterogeneity can be seen as normalizing excitability across neurons. Such network changes echo those resulting from a transition between class I and class II excitability (76), in line with previous experimental observations comparing healthy and epileptogenic neurons (37) as well as other pathologies (77). Such linearization in mean firing rate response is known to have profound consequences on the stability or neural networks (37, 78, 79).

To better understand how this linearization arises, we next evaluated how neuron-to-neuron differences in excitability profiles impacted firing rate covariability. To do so, we considered regimes of stability (in both the homogeneous and heterogeneous cases) and quantified pairwise firing rate correlations between neurons jointly driven by a common modulatory input of varying amplitude. As shown in Fig. 2D, neuron-to-neuron differences in excitability threshold limit covariation between their individual firing rates. We validated these results analytically in the weak connectivity limit (Materials and Methods). These results provide context to those presented in Fig. 1 A and B, in which the steepness of the response of heterogeneous networks is suppressed in favor of smooth uncorrelated fluctuations. The effect was also found to scale with modulatory amplitude. We note that this is a direct consequence of the effect portrayed in Fig. 2C. This analysis suggests that excitability heterogeneity also manifests through the decorrelation of neural responses—thereby limiting network-wide (i.e., correlated) firing rate fluctuations caused by modulatory inputs (59, 80). This decorrelation is complementary yet distinct from the synaptic mechanisms of excitatory–inhibitory balance (81), relying instead on cell-to-cell differences in excitability profiles.

The influence of heterogeneity on stability (and hence volatility) was also confirmed by numerically computing Lyapunov exponents across independent realizations of the network connectivity and independent trials, in which the system exhibits different configurations and is exposed to variable modulatory inputs. The network is unstable whenever Lyapunov exponents are positive, and stable otherwise. As shown in Fig. 2E, persistent positive mean Lyapunov exponents with large variance are observed in the homogeneous case. Such volatility is expected from the theory of nonlinear balanced networks (71, 72). In the heterogeneous case, Lyapunov exponents remained bounded below zero throughout. These results confirm persistent stability and suggest enhanced resilience in the presence of excitability heterogeneity.

Figs. 1 and 2 illustrate the cell and circuit manifestations of the profound dynamical changes that can arise from excitability heterogeneity and its attendant loss. To better understand the underlying mechanisms underlying the cell and circuit effects of excitability heterogeneity on dynamical stability, we harnessed spectral theory for random matrices (28, 61, 82), which is a way of characterizing the stability of large-scale networks and dynamical systems. By construction, our network model is subject to the circular law of random matrix theory (61, 82) in which the eigenvalues—defining stability—are constrained with high probability in a disk in the complex plane. The eigenvalues populating that disk are complex numbers. If they all possess negative real parts, the network is stable: Fluctuations are absorbed by the network, and neural activity returns to its equilibrium after being perturbed. In our network model, such stability manifests through low-intensity, asynchronous neuronal firing smoothly driven by modulatory inputs (e.g., Fig. 2B). However, if some eigenvalues have positive real parts, the network is instead unstable: Neural activity diverges once perturbed, exhibiting synchronous and/or chaotic dynamics (e.g., Fig. 2A). In the intermediate case, when eigenvalues lie close to the imaginary axis, the network is said to reside at a critical point sitting between stability and instability, commonly referred to as metastable. The circular law provides a powerful way to analyze the behavior of a large number of eigenvalues simultaneously, without the challenging task of considering them individually. We direct the reader toward ref. (83) for an excellent introduction and discussion on spectral analysis, the circular law, and its applications in neuroscience.

The circular law (61, 82) has been extensively used in neuroscience to characterize the dynamics of neural circuits (83–86) and in macroecology to examine the stability of complex natural systems, such as food webs (12, 14, 23–28). It provides a powerful and intuitive framework for characterizing resilience in our network model. Given that eigenvalues are constrained within a disk of radius Γ (called spectral radius), one can thus quantify stability by considering how Γ changes with modulatory input, heterogeneity, and other control parameters. The spectral radius may be viewed as a measure of robustness of neuronal populations, indexing how much networks resist change and retain stability.

Geometrically, changes in the spectral radius Γ (reflecting changes in the size of the spectral disk) result either in the clustering or dispersion of eigenvalues centered around the local relaxation gain ℓ; Materials and Methods. As can be seen in Fig. 3A, whenever the spectral radius Γ becomes larger (resp. smaller) than |ℓ|, eigenvalues will cross the imaginary axis, and the network becomes unstable (resp. stable) with high probability (30, 31, 61, 87–89). If the spectral radius remains commensurate with |ℓ|, then the network is considered metastable and in the vicinity of a critical point.

While the net size of the spectral radius determines the system’s stability, how this spectral radius changes with respect to a control parameter (e.g., modulatory input amplitude S_o) reflects the system’s resilience or volatility. If eigenvalues (which are contained in the spectral disk) move significantly in the complex plane due to perturbations (e.g., modulatory input or other control parameters), the network is prone to sudden qualitative changes in dynamics. Indeed, as neural systems can reside in both stable (relaxation) and/or unstable (oscillations, synchrony, chaos) functionally meaningful dynamic regimes, the value of eigenvalues (inferred from the spectral radius Γ) evaluated at a given moment in time conveys little information. It is instead how these eigenvalues change that reflects resilience or volatility. To examine this sensitivity, we subjected the homogeneous network to a thorough spectral analysis (Materials and Methods). By virtue of having homogenous excitability, individual neurons’ steady states were found to be identical across the network and entirely dependent on the modulatory input amplitude (i.e., u_j^o = B + S_o), as expected. This is fully consistent with the dynamics observed in Fig. 2. Over short time scales, the modulatory input S(t)≈S_o can be considered constant: Its influence on the spectral radius Γ may thus be quantified. Indeed, as can be seen in Fig. 3A, both numerically and analytically, the spectral radius was found to be highly sensitive to modulatory input: Changes in S_o resulted in high-amplitude clustering and/or dispersion of the eigenvalues around the relaxation gain, causing frequent transitions between stability and instability. The spectral radius Γ was found to increase with the modulatory input amplitude (S_o), indicating that such fluctuations generally lead to instability.

We confirmed this volatility in Fig. 3B, alongside the alignment between our numerical and analytical calculations. Time-dependent changes in the amplitude of the modulatory input (such as those exemplified in Fig. 2) significantly contract and/or expand spectral radius Γ, whose value intermittently crosses the stability threshold, leading to an alternation between stability and instability. As S_o fluctuates, the network undergoes epochs of instability, alternating with periods where neural activity is either suppressed (S_o strongly inhibiting) and/or saturated (S_o strongly exciting). We note that fast changes in S(t) might cause the network to cross the unstable regime briefly; instability is then difficult to observe since the system does not evolve sufficiently fast to exhibit unstable observable dynamics. Results plotted in Fig. 3C show a high dependence of the spectral radius on modulatory input amplitude (S_o). Stability (i.e., relaxing neural activity, small Γ) characterizes inhibitory and/or low-amplitude modulatory input, while higher amplitudes lead to instability (i.e., divergent, chaotic, and/or synchronous neural activity, large Γ) and eventually saturation (i.e., neural activity plateaus, small Γ). As shown in Fig. 3D, our analysis also revealed that the spectral radius Γ—and hence the dispersion of eigenvalues in the complex plane—increases with network size (N), connection probability (ρ), firing rate response gain (β), as well as net synaptic strength (μ); individually or collectively, all these network features diminish the system’s resilience. This is in line with previous results (25, 30) notably on balanced networks (71, 72), highlighting that homogeneous networks are generically prone to instability. Taken together, our analysis indicates that in sparse balanced and homogeneous networks, high sensitivity to modulatory input and other control parameters underlies the system’s volatility.

Numerous previous studies (39, 41, 60, 75) have shown that heterogeneous neural systems adapt and converge toward a regime of metastability to optimize responses and coding properties. Such metastability manifests itself through critical-like neural activity (90, 91) and/or dynamics residing in the vicinity of a state transition (92). From the perspective of the aforementioned circular law, such dynamical properties emerge whenever these networks are brought toward and operate in dynamical regimes resulting from an intermediate spectral radius, neither too small (i.e., strong stability leading to quiescence; eigenvalues all have negative real parts) nor too large (i.e., strong instability leading to divergence, chaos, and/or synchrony; some eigenvalues have positive real parts).

Our previous findings (37) suggest that intrinsic excitability heterogeneity should improve network resilience, as portrayed in Figs. 1 and 2. To better understand the mechanism behind these dynamics, we adapted the spectral theory for large-scale random systems (28, 61, 82) to expose the influence of excitability heterogeneity on the distribution of eigenvalues. We specifically explored the susceptibility of the spectral radius Γ—and hence the dispersion of eigenvalues in the complex plane—to modulation across various degrees of heterogeneity (σ_H² > 0) (Materials and Methods). Our analysis revealed two main roles played by diversity on network dynamics: a) homeostatic control on network stability and b) the promotion of its resilience.

Indeed, we found that intrinsic excitability heterogeneity is a homeostatic mechanism exerting bidirectional and context-dependent control on network stability: enriching the dynamics whenever they are too poor or, conversely, stabilizing network activity whenever it is too unstable. Indeed, as shown in Fig. 4A, heterogeneity increased the spectral radius (Γ) for small modulatory input amplitudes (S_o). This is also illustrated at the cell and circuit level in Fig. 2B. Excitability heterogeneity maintains network activity during low-drive (S_o) states, and restrains activity in high-drive states. By doing so, heterogeneity prevents both quiescence and excessive firing, while decorrelating firing rates. For low amplitudes of modulation, lack of heterogeneity yields highly stable neural activity that invariably relaxes back to equilibrium whenever perturbed: The spectral radius is infinitesimal, and eigenvalues are clustered around the relaxation gain ℓ. Introducing excitability heterogeneity expands the spectral disk, enriching network dynamics toward instability. Surprisingly, heterogeneity has the opposite effect for higher modulatory input amplitudes, for which the system is highly unstable: Here, heterogeneity instead contracts the spectral disk and stabilizes network dynamics. This contextual control of excitability heterogeneity on stability which depends on modulatory fluctuations suggests that heterogeneity tunes the spectral disk—and hence eigenvalue dispersion—toward an optimal intermediate size.

Indeed, we found that intrinsic excitability heterogeneity is a homeostatic mechanism exerting bidirectional and context-dependent control on network stability: enriching the dynamics whenever they are too poor, or conversely, stabilizing network activity whenever it is too unstable. Indeed, as shown in Fig. 4A, heterogeneity increased the spectral radius (Γ) for small modulatory input amplitudes (S_o). This is also illustrated at the cell and circuit level in Fig. 2B, where excitability heterogeneity maintains network activity during low-drive (S_o) states, while dampening the activity in high-drive states, preventing both quiescence and excessive firing, while decorrelating firing rates. For low amplitudes of modulation, lack of heterogeneity yields highly stable neural activity that invariably relaxes back to equilibrium whenever perturbed: The spectral radius is infinitesimal, and eigenvalues are clustered around the relaxation gain ℓ. Introducing excitability heterogeneity expands the spectral disk, enriching network dynamics toward instability. Surprisingly, higher modulatory input amplitudes, for which the system is highly unstable, led to the opposite. Indeed, heterogeneity was found to here instead contract the spectral disk and stabilize the dynamics. This contextual control of excitability heterogeneity on stability which depends on modulatory fluctuations suggests that heterogeneity tunes the spectral disk—and hence eigenvalue dispersion—toward an optimal intermediate size.

The contextual influence of excitability heterogeneity on network stability described above stems from a damping of spectral radius sensitivity with respect to modulatory input (Fig. 4A). Indeed, sharp changes in Γ caused by S_o (such as those seen in Fig. 3 B and C) were evened out by heterogeneity, resulting in an enrichment or stabilization of the dynamics as the spectral radius is increased or decreased, respectively. Specifically, heterogeneity increased the spectral radius for low modulatory input, while doing the opposite for high amplitudes. The homeostatic influence of heterogeneity on network stability could be confirmed in Fig. 4B. Irrespective of modulatory input amplitude S_o, heterogeneity was found to tune the spectral radius—through either enrichment or stabilization—toward the same intermediate radius.

To confirm the alignment of our mathematical analysis and numerical simulations, we computed the variance of the neuron’s fixed point distribution (i.e., σuo2), which was also found to depend on the degree of heterogeneity (Fig. 4C). Introducing heterogeneity consistently prevented stability transitions, rendering the system more resilient. Indeed, as can be seen Fig. 4D, the transition rate—corresponding to the number of bifurcations observed in the network per unit time—decreased monotonically, confirming the trend seen in Figs. 1 and 2.

Another important conclusion stemming from our analysis is that excitability heterogeneity generically enhances network resilience. As can be seen from Fig. 4, increasing excitability heterogeneity significantly damped spectral radius changes resulting from modulatory input. This implies that excitability heterogeneity anchors eigenvalues in the complex plane, limiting how much they change in response to modulatory input. As discussed before, this limits the probability of changes in stability, preventing bifurcations and confirming resilience. Indeed, excitability heterogeneity made Γ less sensitive to changes in S_o, and by doing so, quenched volatility. This was confirmed in Fig. 5A by systematically varying modulatory input amplitude and the degree of heterogeneity while measuring the spectral radius. We found that heterogeneity damped the sensitivity of the network stability on S_o, as the spectral radius gradually becomes effectively independent of S_o beyond a given degree of heterogeneity (dashed box in Fig. 5A). To encapsulate the effect of excitability heterogeneity on the network’s resilience, we computed both the spectral volatility (κ)—which measures the effective sensitivity of the spectral radius on a given control parameter—as well as the resilience parameter (ℛ)—which is the reciprocal of the spectral volatility—as a function of modulatory input amplitude (i.e., S_o). These metrics quantify how invariant to changes in a given control parameter the eigenvalue distribution is. This is done by looking at variations of the spectral radius, see Materials and Methods. As shown in Fig. 5B, heterogeneity optimized resilience to modulatory input, and the spectral volatility decreased. Collectively, these results demonstrate that excitability heterogeneity greatly enhances the resilience of sparse balanced networks by anchoring the eigenvalues in the complex plane and decreasing the sensitivity of their distribution to modulatory input.

Our results so far indicate that excitability heterogeneity implements a long time-scale homeostatic control mechanism that promotes resilience in networks exposed to modulatory inputs. However, other control parameters might influence the neural systems’ stability over these time scales. Neural systems are subjected to perpetual change, even in the absence of modulatory fluctuations and/or stimuli. Synaptic plasticity is a salient example: During learning, the number of synapses and/or the effective synaptic weights change, as a consequence of processes such as long-term potentiation (LTP) and depression (LTD) (10). Networks undergoing such plasticity-induced structural modifications of their connectivity tend to be weakly resilient. Indeed, most forms of synaptic plasticity lead to the development of instability, in which run-away neural activity departs from baseline and needs to be compensated/stabilized through various homeostatic feedback processes (5, 66, 67), a few of which have found experimental support (93).

We asked whether excitability heterogeneity, on its own, could prevent stability transitions in neural systems undergoing plasticity-induced changes in connectivity. Our previous analysis demonstrates that the spectral radius Γ—and hence the dispersion of eigenvalues in the complex plane—increases with connection probability (ρ) as well as net synaptic strength (μ). This suggests that long time-scale changes in these network features—prone to increase together or independently during learning—generically promote instability and volatility. As can be seen in Fig. 6, this is confirmed numerically: Increasing both the connection probability (Fig. 6A) and synaptic strength (Fig. 6B) over a physiologically realistic range resulted in instability, as measured with the mean Lyapunov exponent. However, this occurred only in the homogeneous case. Indeed, increasing the heterogeneity suppressed this instability with the mean Lyapunov exponent remaining negative over the range of values of explored, i.e., the system did not experience any stability transitions.

Mathematical analysis affirms this finding. Fig. 6C shows that the spectral radius Γ increases monotonically with connection probability (ρ) and synaptic strength (μ) interchangeably, underlying such systems’ volatility and associated vulnerability to stability transitions. Slow time-scale changes in connectivity resulting from plasticity (illustrated by the gray curves linking points a and b in Fig. 6C) result in stability transitions. Introducing heterogeneity moved the effective stability threshold (i.e., Γ(ρ, μ)=|ℓ|) further in parameter space, resulting in overall compensation for the destabilizing influence of increases in connection probability and synaptic strength (c.f., Fig. 3D). In this case, slow time-scale changes in connectivity cause stability transitions to become increasingly unlikely as the net size of the stability region increases. In addition to this stabilizing influence, heterogeneity was also found to promote resilience via anchoring the eigenvalue distribution in the complex plane. We thus computed the resilience metric, now as a function of connection probability (ρ; ℛ_ρ) and synaptic strength (μ; ℛ_μ). As shown in Fig. 6D, increasing the degree of excitability heterogeneity enhanced resilience for both these control parameters, i.e., promoting the persistence of stability by decreasing the spectral volatility and the susceptibility of the spectral radius on changes in connection probability (ρ) and synaptic strength (μ).

We consider a large network of N neurons whose activity evolves according to the interplay between local relaxation, recurrent synaptic connectivity, and slowly varying modulatory input. This model provides a description of dynamics unfolding over extended time scales, hence quantifying mean neuronal activity. The mean somatic membrane potential of neurons u_i(t), i ∈ [1, N] obeys the following set of nonlinear differential equations:

where τ ≫ 1 represents the slow time scale at which the dynamics occur. In the following, we rescale time by t → tτ for convenience. The model in Eq. 1 is both flexible and general, encompassing the mean behavior of a wide scope of interconnected neuron models involving excitatory and inhibitory interactions, such as the celebrated Wilson–Cowan and Jansen–Rit models, for instance. Depending on the spatial scale considered, which remains here undefined, such models can be either considered to be neuron-based (where nodes represent individual neurons—the perspective we adopt here) or populations (where nodes represent assemblies of such neurons), geared toward the characterization of neuronal mean activity across extended time scales. The term $\mathbf{L}[\mathbf{u}]=\ell \mathbf{u}$ is a linear local relaxation term with rate ℓ < 0, and ${(\mathbf{f}[\mathbf{u}])}_i=f(u_i)=\frac{1}{2}(1+\text{erf}[\beta u_i])$ represents the firing rate response function of neurons with the Gaussian error function erf[⋅]. This function relates the membrane potential activity to the firing rate of the neurons. The choice of such response functions is motivated by the fact that they can be generalized to various classes of excitability by a proper adjustment of the response gain β. For instance, small values of β yield class I excitability, while larger values correspond to class II (76). These different cases are both examined in the subsequent analysis. The vector-valued term $\mathbf{H}$ implements node diversity through spatially heterogeneous neuronal excitability, i.e., variable firing rate thresholds between neurons. The entries of $\mathbf{H}$ are sampled from a zero-mean Gaussian distribution of variance ${\sigma }_{\mathbf{H}}^2$. The neuron’s baseline activity $\mathbf{B}<\mathbf{0}$ is a constant offset term used to set the neurons in a subthreshold regime in the absence of input. Last, the network in Eq. 1 is further subjected to a slow modulatory input $\mathbf{S}(t)$.

The connectivity matrix $\mathbf{W}$ in Eq. 1 specifies synaptic coupling between any pair of neurons. We assume randomly distributed excitatory and inhibitory coupling (89), with connection probability ρ. This connectivity motif corresponds to a weighted Erdős–Rényi random graph; we emphasize, however, that the following results may be easily extended to other topologies (e.g., ref. (32)). The strength of these synaptic connections is Gaussian-distributed with mean μ_e and μ_i, variance σ_e² and σ_i², and with probability density functions p_e and p_i, respectively. We ensure that there are no self-connections, i.e., W_ii = 0 ∀i. In addition, we parameterize the relative density of excitatory versus inhibitory connections by a coefficient f, 0 ≤ f ≤ 1. Consequently, the probability density function of synaptic weights W_ij may be written as p = ρfp_e − ρ(1 − f)p_i (89). Moreover, we choose balanced connectivity with $\sum _{j=1}^N{W}_{\mathit{ij}}=0$, i.e., the sum over excitatory and inhibitory synaptic connections vanishes at each node yielding μ_W = 0 and a corresponding synaptic connectivity variance σ_W² (see SI Appendix for details). The mean network activity ⟨u⟩(t) is defined as the average activity across all neurons $⟨u⟩(t)=\sum _{i=1}^N{u}_i(t)/N$.

By construction, the network in Eq. 1 subscribes to the circular law of random matrix theory (61, 82). According to this law, the distribution of eigenvalues—reflecting stability—is constrained with high probability within a disk centered around the local relaxation gain ℓ (called the spectral disk) in the complex plane, whose radius Γ (called the spectral radius) can be determined analytically. If the disk is bounded in the left-hand side of the imaginary axis (i.e., all real parts of these eigenvalues are negative), the network is said to be stable and its activity invariably relaxes back to its equilibrium after a perturbation. In our network model, such stable equilibrium is characterized by weak, asynchronous neuronal firing. If some eigenvalues cross the imaginary axis (i.e., the spectral disk is too large, and some eigenvalues possess positive real parts), the network is said to be unstable, leading to activity that diverges, is synchronous and/or chaotic. In the intermediate case, when dominant eigenvalues (those possessing the largest real part) lie close to the imaginary axis, the network is said to reside at a critical point sitting between stability and instability, commonly referred to as metastable. Through this framework, equilibria of Eq. 1 are stable whenever Γ < |ℓ|, metastable if Γ ≈ |ℓ|, and unstable otherwise. One may hence determine the spectral radius Γ and determine whether stability persists to changes in modulatory input and other control parameters. At short time scales (τ ≪ 1), the slowly varying modulatory input in Eq. 1 can be considered constant, i.e., S(t)=S_o(1, …, 1)^t. Consequently, if all neurons share the same excitability profile, i.e., $\mathbf{H}=0$, we find the stability condition

Eq. 2 further reveals that increasing the network size (N), the connection probability (ρ), the variance of synaptic weights ( σ_W²; implying increases of the mean (μ) and variances σ_e², σ_i²) as well as the firing rate response gain (β) all cause an expansion of the spectral radius Γ, as does the modulatory input amplitude S_o. The dependence of Eq. 2 on these various parameters is plotted in Fig. 3. Consequently, they all lead to instability of Eq. 1, in line with previous studies (32, 110).

The influence of excitability heterogeneity on stability may also be exposed by investigating how diversity in excitability thresholds, i.e., $\mathbf{H}\ne \mathbf{0}$, impacts the spectral radius Γ. Recall that elements of $\mathbf{H}$ are random and sampled from a Gaussian distribution of mean 0 and variance ${\sigma }_{\mathbf{H}}^2$. In this case, the stability condition becomes instead

where $\gamma =1+4{\beta }^2({\sigma }_{{\mathbf{u}}^o}^2+{\sigma }_{\mathbf{H}}^2)$ (see SI Appendix for a detailed derivation).

Resilience refers to the qualitative invariance of dynamical states and the absence of stability transitions. To measure resilience, one may quantify the sensitivity of the spectral radius, i.e., how much Γ fluctuates when exposed to changes to a given control parameter P. To do so, we defined the spectral volatility measure κ_P, reflecting the sensitivity of the spectral disk to change in the parameter P over the range of values that this parameter can take (cf., SI Appendix for more details). Small volatility reflects persistence of the eigenvalue distributions and its overall resistance toward stability transitions. We may thus introduce the reciprocal of κ_P to quantify resilience, i.e., ℛ_P = 1/(1 + κ_P). Note that low (high) volatility corresponds to high (low) resilience.

Given that we have derived the spectral radius analytically in Eq. 3, we are in a position to compute the spectral volatility and resilience analytically. Specifically, if one considers P = S_o, one obtains

where $\gamma =1+4{\beta }^2({\sigma }_{{\mathbf{u}}^o}^2+{\sigma }_{\mathbf{H}}^2)$. Eqs. 4 and 5 show that the volatility and resilience of the network with respect to the modulatory input amplitude both depend on excitability heterogeneity through the factor γ. Whenever σ_H² > 0 increases, the spectral volatility κ_So decreases and the resilience ℛ_So increases.

To generate representative spiking activity in Fig. 2 A and B associated with the dynamics of our network model (i.e., Eq. 1), we simulated Poisson spike trains whose firing rates match those generated by neurons in our network. Specifically, we simulated Poisson processes X_i → Poisson(r_i), where r_i = f_i = f[u_i + h_i] corresponds to the firing rate of individual neurons. The resulting exemplar spiking activity corresponds to a given a realization of such Poisson processes whose rates are identical to the ones predicted by our model.

Normalized mean network firing rates $\overline{r}$ in Fig. 2C can be computed with and without heterogeneity. Since $\overline{r}=\sum _{i=1}^N{f}_i/N$, in the N ≫ 1 limit of large networks, we find

see SI Appendix for more details.

Excitability heterogeneity also shapes pairwise firing rate correlations. Over long time scales (i.e., τ ≫ 1) the modulatory input is not constant and may be considered a normally distributed and uncorrelated random process of mean 0 and variance 2S_o², resulting in correlated fluctuations in neuronal activity. One may further consider the weak coupling limit where the recurrent connectivity is small compared to modulatory input amplitude (σ_W² ≪ S_o²). In this regime, the activity of individual neurons in Eq. 1 is strongly driven by modulatory inputs, exhibiting a normally distributed probability density function of mean 0 and variance S_o². We note that these approximations hold whenever the network is stable (i.e., Γ < |ℓ|).

For neurons i and j with excitability thresholds h_i and h_j, respectively, the mean pairwise correlation coefficient c_ij between the firing rates r_i = f[u_i + h_i] and r_j = f[u_j + h_j] may be expressed as

where cov[r_ir_j](Δh_ij) corresponds to the covariance expressed as a function of the difference Δh_ij = h_j − h_i in excitability between two neurons i and j, and Var[r_{i, j}] is the firing rate variance (see SI Appendix for further details).

A.H., S.R., T.A.V., and J.L. designed research; A.H. and J.L. performed research; A.H. and J.L. analyzed data; and A.H., S.R., T.A.V., and J.L. wrote the paper.

The authors declare no competing interest.