<figure><img src="image_1.png" /></figure>

HAL
open science

# Diversity-induced trivialization and resilience of neural dynamics

Axel Hutt, Daniel Trotter, Aref Pariz, Taufik A Valiante, Jérémie Lefebvre

► To cite this version:

Axel Hutt, Daniel Trotter, Aref Pariz, Taufik A Valiante, Jérémie Lefebvre. Diversity-induced trivialization and resilience of neural dynamics. Chaos: An Interdisciplinary Journal of Nonlinear Science, In press. hal-04162586v2

HAL Id: hal-04162586

https://inria.hal.science/hal-04162586v2

Submitted on 23 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Diversity-induced trivialization and resilience of neural dynamics

Axel Hutt, Daniel Trotter, Aref Pariz, Taufik A Valiante, Jérémie Lefebvre

► To cite this version:

Axel Hutt, Daniel Trotter, Aref Pariz, Taufik A Valiante, Jérémie Lefebvre. Diversity-induced trivialization and resilience of neural dynamics. Chaos: An Interdisciplinary Journal of Nonlinear Science, 2024, 34 (1), 10.1063/5.0165773. hal-04162586v2

HAL Id: hal-04162586

https://inria.hal.science/hal-04162586v2

Submitted on 23 Feb 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Diversity-induced trivialization and resilience of neural dynamics

Axel Hutt,<fnref n="1" /> Daniel Trotter,<fnref n="2" />,<fnref n="3" /> Aref Pariz,<fnref n="3" />,<fnref n="4" /> Taufik A Valiante,<fnref n="3" />,<fnref n="5" /> and Jérémie Lefebvre<fnref n="2" />, <fnref n="3" />, <fnref n="4" />, <fnref n="6" />

Heterogeneity is omnipresent in living systems and biophysical diversity enriches the systems' dynamical repertoire. However, it remains challenging to reconcile with the robustness and persistence of the systems functions over time, which is called resilience. To better understand the underlying mechanism of resilience, we considered a nonlinear neural network model focussing on the relationship between excitability heterogeneity of neurons and resilience. To quantify the degree of resilience, we considered the number of stationary states present in the system and how they are affected by parameters. This impact on the number of stationary states is known as trivialization. We analyzed both analytically and numerically gradient and non-gradient systems modeled as non-linear sparse neural networks evolving over long time scales. Excitability heterogeneity in neurons tuned network stability in a context-dependent way, quenched the number of stationary states and enhanced resilience. This heterogeneity was found to implement a homeostatic control mechanism enhancing network resilience to changes in network size and connection probability by quenching the system's dynamic volatility.

Recent research has identified large diversity in living systems, such as variations in bio-physical features in the brain. This raises the question on the functional role of such a diversity and how it impacts the systems stability. By mathematical analysis and computational simulations of trivialization in a sparse random neural network, we demonstrate that heterogeneity in neuron excitability represents a homeostatic control mechanism. Hence excitability heterogeneity tunes the system's resilience to perturbations. This highlights the importance of bio-physical diversity in the persistence of a system's function.

## I. INTRODUCTION

Complex biological systems, such as the neural circuits of the brain, are surprisingly resilient over time<fnref n="1" />. The dynamics of these networks resist change, persisting in the face of perturbations, such as neuromodulation, stimuli, pathological insults, and plasticity occurring during development and learning. Failure to compensate for such perturbations would otherwise undermine the reliability of brain function and predispose these circuits to pathological states, such as those seen in epilepsy<fnref n="2" />. What makes those systems so stable in spite of their complexity? Identifying the mechanisms and/or principles responsible for such robustness is a formidable experimental and math-

ematical problem with important implications for brain function and disease. Enduring stability is a desirable, yet challenging trait to implement in engineering. Optimizing organizational and societal designs for efficiency often jeopardizes resilience, leading to fragile systems prone to collapse.<fnref n="3" />⁻<fnref n="6" />

Despite important advances in our understanding of complex systems stability<fnref n="3" />,<fnref n="5" />,⁷⁻¹⁷, the manifest resilience of neural circuits remains difficult to reconcile with their widespread and constantly evolving diversity. Efforts in identifying the features imparting resilience to neural networks have historically been oriented towards their topology (i.e. the properties of their connectivity), yet results remain perplexing⁸⁻¹⁰,¹⁶,¹⁸. Indeed, the hallmark of complex systems is their high susceptibility to instability: traits that accompany complexity, such as increased dimensionality⁹,¹⁸⁻²⁰, strength of connections¹⁰⁻¹²,¹⁴,¹⁸,²¹, motifs richness¹⁰⁻¹²,¹⁴,¹⁸,²¹, connectance⁸,⁹,¹⁴,¹⁸,²¹ and degree heterogeneity²² generally make complex systems fragile<fnref n="4" /> and prone to instability through a so-called May-Wigner state transition⁸⁻¹⁰,¹⁴,¹⁵,²³⁻²⁵. Indeed, increases in such topological characteristics most often hinder resilience instead of promoting it<fnref n="5" />,⁷,¹⁷,¹⁸,²⁶,²⁷. In light of the constant reorganization and complexification during development and learning, neural circuits require compensatory strategies to preserve the robustness of brain dynamics and function over time.

One way to quantify the resilience of a given dynamical system is to measure the number of its equilibria²³,²⁸⁻³⁰.

<footnote n="1">¹) MLMS, MIMESIS, Université de Strasbourg, CNRS, Inria, ICube, 67000 Strasbourg,
France</footnote>

<footnote n="2">²) Department of Physics, University of Ottawa, Ottawa, ON K1N 6N5, Canada</footnote>

<footnote n="3">³) Krembil Brain Institute, University Health Network, Toronto, ON M5T 0S8,
Canada</footnote>

<footnote n="4">⁴) Department of Biology, University of Ottawa, Ottawa, ON K1N 6N5, Canada</footnote>

<footnote n="5">⁵) Department of Electrical and Computer Engineering; Institute of Medical Science;
Institute of Biomedical Engineering; University of Toronto, ;Division of Neurosurgery,
Department of Surgery; CRANIA (Center for Advancing Neurotechnological Innovation to Application;
Max Planck-University of Toronto Center for Neural Science and Technology, Toronto, ON,
Canada</footnote>

<footnote n="6">⁶) Department of Mathematics, University of Toronto, Toronto, ON M5S 2E4, Canada</footnote>
# Diversity-induced trivialization and resilience of neural dynamics

Axel Hutt,¹ Daniel Trotter,²,³ Aref Pariz,³,⁴ Taufik A Valiante,³,⁵ and Jérémie Lefebvre², ³, ⁴, ⁶

¹) MLMS, MIMESIS, Université de Strasbourg, CNRS, Inria, ICube, 67000 Strasbourg,
France

²) Department of Physics, University of Ottawa, Ottawa, ON K1N 6N5, Canada

³) Krembil Brain Institute, University Health Network, Toronto, ON M5T 0S8,
Canada

⁴) Department of Biology, University of Ottawa, Ottawa, ON K1N 6N5, Canada

⁵) Department of Electrical and Computer Engineering; Institute of Medical Science;
Institute of Biomedical Engineering; University of Toronto, ;Division of Neurosurgery,
Department of Surgery; CRANIA (Center for Advancing Neurotechnological Innovation to Application;
Max Planck-University of Toronto Center for Neural Science and Technology, Toronto, ON,
Canada

⁶) Department of Mathematics, University of Toronto, Toronto, ON M5S 2E4, Canada

Heterogeneity is omnipresent in living systems and biophysical diversity enriches the systems' dynamical repertoire. However, it remains challenging to reconcile with the robustness and persistence of the systems functions over time, which is called resilience. To better understand the underlying mechanism of resilience, we considered a nonlinear neural network model focussing on the relationship between excitability heterogeneity of neurons and resilience. To quantify the degree of resilience, we considered the number of stationary states present in the system and how they are affected by parameters. This impact on the number of stationary states is known as trivialization. We analyzed both analytically and numerically gradient and non-gradient systems modeled as non-linear sparse neural networks evolving over long time scales. Excitability heterogeneity in neurons tuned network stability in a context-dependent way, quenched the number of stationary states and enhanced resilience. This heterogeneity was found to implement a homeostatic control mechanism enhancing network resilience to changes in network size and connection probability by quenching the system's dynamic volatility.

> Recent research has identified large diversity in living systems, such as variations in bio-physical features in the brain. This raises the question on the functional role of such a diversity and how it impacts the systems stability. By mathematical analysis and computational simulations of trivialization in a sparse random neural network, we demonstrate that heterogeneity in neuron excitability represents a homeostatic control mechanism. Hence excitability heterogeneity tunes the system's resilience to perturbations. This highlights the importance of bio-physical diversity in the persistence of a system's function.

## I. INTRODUCTION

Complex biological systems, such as the neural circuits of the brain, are surprisingly resilient over time¹. The dynamics of these networks resist change, persisting in the face of perturbations, such as neuromodulation, stimuli, pathological insults, and plasticity occurring during development and learning. Failure to compensate for such perturbations would otherwise undermine the reliability of brain function and predispose these circuits to pathological states, such as those seen in epilepsy². What makes those systems so stable in spite of their complexity? Identifying the mechanisms and/or principles responsible for such robustness is a formidable experimental and math-

ematical problem with important implications for brain function and disease. Enduring stability is a desirable, yet challenging trait to implement in engineering. Optimizing organizational and societal designs for efficiency often jeopardizes resilience, leading to fragile systems prone to collapse.³⁻⁶

Despite important advances in our understanding of complex systems stability³,⁵,⁷⁻¹⁷, the manifest resilience of neural circuits remains difficult to reconcile with their widespread and constantly evolving diversity. Efforts in identifying the features imparting resilience to neural networks have historically been oriented towards their topology (i.e. the properties of their connectivity), yet results remain perplexing⁸⁻¹⁰,¹⁶,¹⁸. Indeed, the hallmark of complex systems is their high susceptibility to instability: traits that accompany complexity, such as increased dimensionality⁹,¹⁸⁻²⁰, strength of connections¹⁰⁻¹²,¹⁴,¹⁸,²¹, motifs richness¹⁰⁻¹²,¹⁴,¹⁸,²¹, connectance⁸,⁹,¹⁴,¹⁸,²¹ and degree heterogeneity²² generally make complex systems fragile⁴ and prone to instability through a so-called May-Wigner state transition⁸⁻¹⁰,¹⁴,¹⁵,²³⁻²⁵. Indeed, increases in such topological characteristics most often hinder resilience instead of promoting it⁵,⁷,¹⁷,¹⁸,²⁶,²⁷. In light of the constant reorganization and complexification during development and learning, neural circuits require compensatory strategies to preserve the robustness of brain dynamics and function over time.

One way to quantify the resilience of a given dynamical system is to measure the number of its equilibria²³,²⁸⁻³⁰.
The number of fixed points is indeed directly linked to a system's susceptibility to multistability and potential for bifurcations, reflecting a predisposition to transitions between states in the wake of disturbances. One may indeed interpret the dynamics of complex systems as evolving on a (a *priorily* unknown) high-dimensional landscape, occupied by various optima, corresponding to either stable or unstable fixed points. The topology of this landscape - that is the number of equilibria and the size of their respective basins of attraction - determines the response of such systems to perturbations. Resilience decreases with increasing number of fixed points: the more equilibria are present, the more complex the landscape is, predisposing this system to transitions between qualitatively different dynamics. Seminal studies<sup>23,28</sup> have shown that, in general, neural networks are susceptible to May-Wigner transitions in which the number of equilibria populating such landscapes increases exponentially as they transition towards instability, notably as a result of increase in connectivity<sup>9,10,12,17,27,31</sup>. The existence and ubiquity of such transitions, commonly accompanying the emergence of chaotic dynamics<sup>32</sup>, are *a priori* difficult to reconcile with the manifest resilience and robustness of neural circuits.

For many of these studies, however, one central assumption is that individual neurons, or nodes of these neural networks, are either indistinguishable from each other or alternatively break down into discrete classes (e.g. excitatory pyramidal cells, inhibitory interneurons) sharing similar excitability profiles and responses to stimuli. However, biophysical heterogeneity amongst neurons is one of the most salient features of neural tissue<sup>33-37</sup>, even within the same cell type<sup>38-43</sup>, whose phenotypes are far from identical and instead distributed along a continuum. Pioneering experimental<sup>2,44-52</sup> and theoretical<sup>2,53-61</sup> studies have explored such heterogeneity, especially intrinsic variability in excitability, and demonstrated that it imparts neural circuits with optimized coding capabilities<sup>54,56,58,59,62</sup> and robustness<sup>47,48,52,63,64</sup>, whose dynamics remain qualitatively invariant over a wide range neuromodulatory perturbations and/or disturbances in connectivity<sup>65</sup>. Notably, excitability heterogeneity - a feature originating from various combinations of single cell biophysical attributes - has been shown to stabilize neural dynamics away from pathological brain states<sup>2</sup>.

These results echo the conclusions of numerous studies in the physical<sup>66-69</sup> and biological<sup>70,71</sup> sciences demonstrating that noise and/or disorder may help quench complexity to support reliability<sup>65,72</sup>. Such studies make an eloquent argument as to the counter intuitive role of randomness and variability in promoting the predictability and reliability of complex systems dynamics. Interestingly, statistical mechanics approaches have also convincingly demonstrated that under specific conditions, nonlinear systems may simplify through so called *topological trivialization*<sup>24,25</sup> hereafter called trivialization. This

term, introduced by Fyodorov and Le Doussal, refers to the process by which the total number of equilibria - and hence net complexity of the aforementioned landscape - is quenched. An eloquent example for this phenomenon is the special case of gradient systems, whose dynamics evolve in highly non-linear and frequently random potentials. Through trivialization, the geometry of such potentials is restructured and its quadratic component enhanced, constraining the dynamics around a limited set of equilibria<sup>23,24</sup>. As a corollary, the net stability and resilience of such system are hence both reinforced. Interestingly, trivialization has been shown to result from intrinsic heterogeneity in spin glass systems<sup>29</sup>, suggesting that diversity in complex biological system might serve a similar functional role, limiting the dynamical repertoire of otherwise highly unstable systems to promote their resilience.

Inspired by these results, we here explore this hypothesis from the perspective of heterogeneity, and how it might enhance the resilience of large-scale neural networks to perturbations through trivialization. That is, intrinsic neural diversity may promote stability by reducing the number of fixed points of a generic class of neural networks, thereby counter-intuitively quenching their dynamical repertoire. We first summarize the elegant theoretical framework developed<sup>23,25,28,30</sup> to quantify the number of equilibria of a system susceptible to May-Wigner stability transitions. Using this framework we demonstrate how the dispersion in excitability profile delays such transitions while suppressing the number of equilibria through trivialization. We discuss these results and explore their potential functional implications, notably in stabilizing neural dynamics during development and learning.

## II. TRIVIALIZATION: THEORY

### A. Neural Network Model

Determining the number of fixed points for high-dimensional non-linear systems is a notoriously challenging problem<sup>23-25,28,30</sup>, that possesses a large range of applications in optimization, physics and biology. Fortunately, such calculations remain tractable given some simplifying assumptions. Our analysis is motivated by a widely used type of neural network<sup>32,73</sup> in which the activity $u_j$ of a neuron obeys the following set of coupled evolution equations

$$ \frac{d}{dt} u_i(t) = L[u_i] + \sum_{j=1}^{N} w_{ij} f_j[u_j] \quad (1) $$

where $L[u] = -u$ is a linear relaxation operator and $f_j[u] = \frac{1}{2}(1 + \text{erf}[\beta(u + h_j)])$ is a nonlinear sigmoidal activation function with gain $\beta$, where $h_j$ sets the individual excitability of neurons. The synaptic weights $w_{ij}$
The number of fixed points is indeed directly linked to a system's susceptibility to multistability and potential for bifurcations, reflecting a predisposition to transitions between states in the wake of disturbances. One may indeed interpret the dynamics of complex systems as evolving on a (a *priorily* unknown) high-dimensional landscape, occupied by various optima, corresponding to either stable or unstable fixed points. The topology of this landscape - that is the number of equilibria and the size of their respective basins of attraction - determines the response of such systems to perturbations. Resilience decreases with increasing number of fixed points: the more equilibria are present, the more complex the landscape is, predisposing this system to transitions between qualitatively different dynamics. Seminal studies<sup>23,28</sup> have shown that, in general, neural networks are susceptible to May-Wigner transitions in which the number of equilibria populating such landscapes increases exponentially as they transition towards instability, notably as a result of increase in connectivity<sup>9,10,12,17,27,31</sup>. The existence and ubiquity of such transitions, commonly accompanying the emergence of chaotic dynamics<sup>32</sup>, are *a priori* difficult to reconcile with the manifest resilience and robustness of neural circuits.

For many of these studies, however, one central assumption is that individual neurons, or nodes of these neural networks, are either indistinguishable from each other or alternatively break down into discrete classes (e.g. excitatory pyramidal cells, inhibitory interneurons) sharing similar excitability profiles and responses to stimuli. However, biophysical heterogeneity amongst neurons is one of the most salient features of neural tissue<sup>33-37</sup>, even within the same cell type<sup>38-43</sup>, whose phenotypes are far from identical and instead distributed along a continuum. Pioneering experimental<sup>2,44-52</sup> and theoretical<sup>2,53-61</sup> studies have explored such heterogeneity, especially intrinsic variability in excitability, and demonstrated that it imparts neural circuits with optimized coding capabilities<sup>54,56,58,59,62</sup> and robustness<sup>47,48,52,63,64</sup>, whose dynamics remain qualitatively invariant over a wide range neuromodulatory perturbations and/or disturbances in connectivity<sup>65</sup>. Notably, excitability heterogeneity - a feature originating from various combinations of single cell biophysical attributes - has been shown to stabilize neural dynamics away from pathological brain states<sup>2</sup>.

These results echo the conclusions of numerous studies in the physical<sup>66-69</sup> and biological<sup>70,71</sup> sciences demonstrating that noise and/or disorder may help quench complexity to support reliability<sup>65,72</sup>. Such studies make an eloquent argument as to the counter intuitive role of randomness and variability in promoting the predictability and reliability of complex systems dynamics. Interestingly, statistical mechanics approaches have also convincingly demonstrated that under specific conditions, nonlinear systems may simplify through so called *topological trivialization*<sup>24,25</sup> hereafter called trivialization. This

term, introduced by Fyodorov and Le Doussal, refers to the process by which the total number of equilibria - and hence net complexity of the aforementioned landscape - is quenched. An eloquent example for this phenomenon is the special case of gradient systems, whose dynamics evolve in highly non-linear and frequently random potentials. Through trivialization, the geometry of such potentials is restructured and its quadratic component enhanced, constraining the dynamics around a limited set of equilibria<sup>23,24</sup>. As a corollary, the net stability and resilience of such system are hence both reinforced. Interestingly, trivialization has been shown to result from intrinsic heterogeneity in spin glass systems<sup>29</sup>, suggesting that diversity in complex biological system might serve a similar functional role, limiting the dynamical repertoire of otherwise highly unstable systems to promote their resilience.

Inspired by these results, we here explore this hypothesis from the perspective of heterogeneity, and how it might enhance the resilience of large-scale neural networks to perturbations through trivialization. That is, intrinsic neural diversity may promote stability by reducing the number of fixed points of a generic class of neural networks, thereby counter-intuitively quenching their dynamical repertoire. We first summarize the elegant theoretical framework developed<sup>23,25,28,30</sup> to quantify the number of equilibria of a system susceptible to May-Wigner stability transitions. Using this framework we demonstrate how the dispersion in excitability profile delays such transitions while suppressing the number of equilibria through trivialization. We discuss these results and explore their potential functional implications, notably in stabilizing neural dynamics during development and learning.

## II. TRIVIALIZATION: THEORY

### A. Neural Network Model

Determining the number of fixed points for high-dimensional non-linear systems is a notoriously challenging problem<sup>23-25,28,30</sup>, that possesses a large range of applications in optimization, physics and biology. Fortunately, such calculations remain tractable given some simplifying assumptions. Our analysis is motivated by a widely used type of neural network<sup>32,73</sup> in which the activity $u_j$ of a neuron obeys the following set of coupled evolution equations

$$ \frac{d}{dt} u_i(t) = L[u_i] + \sum_{j=1}^{N} w_{ij} f_j[u_j] \quad (1) $$

where $L[u] = -u$ is a linear relaxation operator and $f_j[u] = \frac{1}{2}(1 + \text{erf}[\beta(u + h_j)])$ is a nonlinear sigmoidal activation function with gain $\beta$, where $h_j$ sets the individual excitability of neurons. The synaptic weights $w_{ij}$
are random, normally distributed with mean $\mu_W$ and variance $\sigma_w^2 > 0$ and connection probability $\rho$. One may further consider balanced connectivity$^{73}$, in which strength of these synaptic connections are individually Gaussian-distributed with mean $\mu_e$ and $\mu_i$, variance $\sigma_e^2$ and $\sigma_i^2$ and with probability density functions $p_e(x)$ and $p_i(x)$, respectively, with $x \in \mathbb{R}$. The relative density of excitatory versus inhibitory connections is parameterized by $g$, $0 \le g \le 1$, so that the probability density function of synaptic weights $w_{ij}$ may be written as

$$p(x) = \rho g p_e(x) - \rho(1-g)p_i(x).$$

Based on these criteria, the mean and variance of the synaptic connectivity reads$^{73}$

$$\mu_w \equiv E[w_{ij}] = \rho(g\mu_e - (1-g)\mu_i)$$

$$\sigma_w^2 \equiv \mathrm{Var}[w_{ij}] = \rho(g\sigma_e^2 + (1-g)\sigma_i^2 + g\mu_w^2/(1-g)).$$

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. see the work of Yann et al.$^{22}$.

The neural network, whose dynamics are described by Eq. (1), has been extensively studied for vanishing average connectivity ($\mu_w = 0$) in the homogeneous case, i.e. neurons with identical excitability profiles ($h_j = 0$). Such systems are well known to exhibit May-Wigner-type transitions$^{32,73–75}$ and an associated exponential increase in the number of equilibria$^{23,28}$ whenever the spectral radius $\Gamma$ becomes larger than unity, i.e. $\Gamma > 1$. Note that a detailed calculation for $\Gamma$ will be provided below. This phenomenon echoes numerous seminal studies linking complexity and stability in high-dimensional biological networks$^{3,5,7–11,14–17,19,21,22,26,27,31}$ whose results hold on the spectral properties of random Gaussian matrices$^{76–78}$.

<figure><img src="image_3.png" /><figcaption>**Figure 1.** Example of trivialization in a gradient system. A The dynamics of the system in Eq. 3 are constrained within a potential V whose shape changes with excitability heterogeneity. In absence of heterogeneity ($\sigma_h^2 = 0$; black line) the potential is bimodal, reflecting the multistability of Eq. 3 and its susceptibility to state transitions. Increase in heterogeneity ($\sigma_h^2 > 0$; blue lines) causes the potential V to become purely quadratic. B Reduction of the number of fixed points of Eq. 3 as excitability heterogeneity is increased. Such fixed points correspond to the (vanishing) maxima of the potential V seen in panel A. Parameters are given $\beta = 15$, $x^o = 0.6$, $\mu_h = 0.25$ and $\sigma_w^2 = 0$.</figcaption></figure>

## B. Gradient systems

Diversity-induced trivialization can be easily observed in simple gradient networks. Let us consider the ensemble average of the system state $\bar{u} \equiv \langle u_i \rangle_N = \sum_{i=1}^N u_i / N$. Then Eq. (1) yields

$$\frac{1}{N} \sum_{i=1}^{N} \frac{du_i}{dt} = \frac{1}{N} \sum_{i=1}^{N} L[u_i] + \frac{1}{N} \sum_{i,j=1}^{N} w_{ij} f_j[u_j]$$

$$\frac{d\bar{u}}{dt} = L[\bar{u}] + \sum_{j=1}^{N} w_j f_j[u_j]$$

with the mean connectivity at each node $W_j = \sum_{i=1}^N w_{ij} / N$. Assuming that firing rates (i.e. $f_j[u_j]$) are independent of synaptic weights (i.e. $w_{ij}$), which holds only whenever $N$ is large$^{65,79,80}$, and by virtue of the network connectivity properties defined above, $W_j = \mu_w + X_j$ with random perturbations $X_j \sim \mathcal{N}(0, \sigma_w^2 / N)$ yields

$$\begin{aligned} \frac{d\bar{u}}{dt} &= L[\bar{u}] + \sum_{j=1}^{N} (\mu_w + X_j) f_j[u_j] \\ &= L[\bar{u}] + N\mu_w \langle f_j \rangle_N + N \langle X_j f_j[u_j] \rangle_N . \end{aligned} \quad (2)$$

Since, in general, $\langle YW \rangle_N = \langle Y \rangle_N \langle W \rangle_N + S$ and $S \sim \mathcal{N}(0, \mathrm{Var}[YW]/N)$ with $\mathrm{Var}[YW] = \mathrm{Var}[Y]\mathrm{Var}[W] + \mathrm{Var}[Y]\langle W \rangle_N^2 + \mathrm{Var}[W]\langle Y \rangle_N^2$ for two statistically independent random variables $Y, W$, Eq. (2) reads

$$\frac{d\bar{u}}{dt} = L[\bar{u}] + N\mu_w \langle f_j \rangle_N + \sigma_w \sqrt{\mathrm{Var}[f_j] + \langle f_j \rangle_N^2} Z$$

with the random variable $Z \sim \mathcal{N}(0, 1)$.

Assuming that the excitability $h_j$ of individual neurons in Eq. (1) are normally distributed with mean $\mu_h$ and variance $\sigma_h^2$, one may consider setting $w_{ij} = x_{ij} / N$ yielding $\mu_w = x^o / N$ and $\sigma_w^2 = \sigma_x^2 / N^2$. Then in the limit $N \to \infty$,$^{14,79–81}$ Eq. (1) admits a well-defined mean field

$$\frac{d}{dt}\bar{u}(t) = L[\bar{u}] + x^o F[\bar{u}]. \quad (3)$$

where $F[\bar{u}]$ corresponds to the statistical average

$$\begin{aligned} \langle f_j \rangle_N &\approx F[\bar{u}] = \int_{-\infty}^{\infty} \frac{1}{2} (1 + \mathrm{erf}(\beta(\bar{u} + h))) p_h(h) dh \\ &= \frac{1}{2} \left( 1 + \mathrm{erf} \left[ \frac{\beta(\bar{u} + \mu_h)}{\gamma} \right] \right) \end{aligned} \quad (4)$$

with $p_h(h)$ as the Gaussian distribution with mean $\mu_h$ and variance $\sigma_h^2$ and $\gamma = \sqrt{1 + 2\beta^2\sigma_h^2}$. Equation (4) holds true whenever $N$ is sufficiently large$^{82,83}$.

Now Eq. (3) may be written in terms of the gradient of a scalar potential $V(\bar{u})$, i.e

$$\frac{d}{dt}\bar{u}(t) = -\frac{dV(\bar{u})}{d\bar{u}}$$
are random, normally distributed with mean $\mu_W$ and variance $\sigma_w^2 > 0$ and connection probability $\rho$. One may further consider balanced connectivity$^{73}$, in which strength of these synaptic connections are individually Gaussian-distributed with mean $\mu_e$ and $\mu_i$, variance $\sigma_e^2$ and $\sigma_i^2$ and with probability density functions $p_e(x)$ and $p_i(x)$, respectively, with $x \in \mathbb{R}$. The relative density of excitatory versus inhibitory connections is parameterized by $g$, $0 \le g \le 1$, so that the probability density function of synaptic weights $w_{ij}$ may be written as

$$p(x) = \rho g p_e(x) - \rho(1-g)p_i(x).$$

Based on these criteria, the mean and variance of the synaptic connectivity reads$^{73}$

$$\mu_w \equiv E[w_{ij}] = \rho(g\mu_e - (1-g)\mu_i)$$

$$\sigma_w^2 \equiv \mathrm{Var}[w_{ij}] = \rho(g\sigma_e^2 + (1-g)\sigma_i^2 + g\mu_w^2/(1-g)).$$

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. see the work of Yann et al.$^{22}$.

The neural network, whose dynamics are described by Eq. (1), has been extensively studied for vanishing average connectivity ($\mu_w = 0$) in the homogeneous case, i.e. neurons with identical excitability profiles ($h_j = 0$). Such systems are well known to exhibit May-Wigner-type transitions$^{32,73–75}$ and an associated exponential increase in the number of equilibria$^{23,28}$ whenever the spectral radius $\Gamma$ becomes larger than unity, i.e. $\Gamma > 1$. Note that a detailed calculation for $\Gamma$ will be provided below. This phenomenon echoes numerous seminal studies linking complexity and stability in high-dimensional biological networks$^{3,5,7–11,14–17,19,21,22,26,27,31}$ whose results hold on the spectral properties of random Gaussian matrices$^{76–78}$.

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Example of trivialization in a gradient system. A The dynamics of the system in Eq. 3 are constrained within a potential $V$ whose shape changes with excitability heterogeneity. In absence of heterogeneity ($\sigma_h^2 = 0$; black line) the potential is bimodal, reflecting the multistability of Eq. 3 and its susceptibility to state transitions. Increase in heterogeneity ($\sigma_h^2 > 0$; blue lines) causes the potential $V$ to become purely quadratic. B Reduction of the number of fixed points of Eq. 3 as excitability heterogeneity is increased. Such fixed points correspond to the (vanishing) maxima of the potential $V$ seen in panel A. Parameters are given $\beta = 15$, $x^o = 0.6$, $\mu_h = 0.25$ and $\sigma_w^2 = 0$.</figcaption></figure>

## B. Gradient systems

Diversity-induced trivialization can be easily observed in simple gradient networks. Let us consider the ensemble average of the system state $\bar{u} \equiv \langle u_i \rangle_N = \sum_{i=1}^N u_i / N$. Then Eq. (1) yields

$$\frac{1}{N} \sum_{i=1}^{N} \frac{du_i}{dt} = \frac{1}{N} \sum_{i=1}^{N} L[u_i] + \frac{1}{N} \sum_{i,j=1}^{N} w_{ij} f_j[u_j] \\ \frac{d\bar{u}}{dt} = L[\bar{u}] + \sum_{j=1}^{N} w_j f_j[u_j]$$

with the mean connectivity at each node $W_j = \sum_{i=1}^N w_{ij} / N$. Assuming that firing rates (i.e. $f_j[u_j]$) are independent of synaptic weights (i.e. $w_{ij}$), which holds only whenever $N$ is large$^{65,79,80}$, and by virtue of the network connectivity properties defined above, $W_j = \mu_w + X_j$ with random perturbations $X_j \sim \mathcal{N}(0, \sigma_w^2 / N)$ yields

$$\begin{aligned} \frac{d\bar{u}}{dt} &= L[\bar{u}] + \sum_{j=1}^{N} (\mu_w + X_j) f_j[u_j] \\ &= L[\bar{u}] + N\mu_w \langle f_j \rangle_N + N \langle X_j f_j[u_j] \rangle_N . \end{aligned} \quad (2)$$

Since, in general, $\langle YW \rangle_N = \langle Y \rangle_N \langle W \rangle_N + S$ and $S \sim \mathcal{N}(0, \mathrm{Var}[YW]/N)$ with $\mathrm{Var}[YW] = \mathrm{Var}[Y]\mathrm{Var}[W] + \mathrm{Var}[Y]\langle W \rangle_N^2 + \mathrm{Var}[W]\langle Y \rangle_N^2$ for two statistically independent random variables $Y, W$, Eq. (2) reads

$$\frac{d\bar{u}}{dt} = L[\bar{u}] + N\mu_w \langle f_j \rangle_N + \sigma_w \sqrt{\mathrm{Var}[f_j] + \langle f_j \rangle_N^2} Z$$

with the random variable $Z \sim \mathcal{N}(0, 1)$.

Assuming that the excitability $h_j$ of individual neurons in Eq. (1) are normally distributed with mean $\mu_h$ and variance $\sigma_h^2$, one may consider setting $w_{ij} = x_{ij} / N$ yielding $\mu_w = x^o / N$ and $\sigma_w^2 = \sigma_x^2 / N^2$. Then in the limit $N \to \infty$,$^{14,79–81}$ Eq. (1) admits a well-defined mean field

$$\frac{d}{dt}\bar{u}(t) = L[\bar{u}] + x^o F[\bar{u}]. \quad (3)$$

where $F[\bar{u}]$ corresponds to the statistical average

$$\begin{aligned} \langle f_j \rangle_N &\approx F[\bar{u}] = \int_{-\infty}^{\infty} \frac{1}{2} (1 + \mathrm{erf}(\beta(\bar{u} + h))) p_h(h) dh \\ &= \frac{1}{2} \left( 1 + \mathrm{erf} \left[ \frac{\beta(\bar{u} + \mu_h)}{\gamma} \right] \right) \end{aligned} \quad (4)$$

with $p_h(h)$ as the Gaussian distribution with mean $\mu_h$ and variance $\sigma_h^2$ and $\gamma = \sqrt{1 + 2\beta^2\sigma_h^2}$. Equation (4) holds true whenever $N$ is sufficiently large$^{82,83}$.

Now Eq. (3) may be written in terms of the gradient of a scalar potential $V(\bar{u})$, i.e

$$\frac{d}{dt}\bar{u}(t) = -\frac{dV(\bar{u})}{d\bar{u}}$$
defined by

$$
V(\bar{u}) = \frac{\bar{u}^2}{2} - \frac{\chi^o}{2\sqrt{\pi}\beta} \left( \gamma e^{\frac{-\beta^2(\mu_h - \bar{u})^2}{\gamma}} \right. \\
\left. + \sqrt{\pi}\beta \left( (\mu_h - \bar{u})\text{erf}\left[ \frac{\beta(\mu_h - \bar{u})}{\sqrt{1 + 2\beta^2\sigma_h^2}} \right] + \bar{u} \right) \right) \quad (5)
$$

Figure 1 depicts trivialization that occurs in Eq. (3) once excitability heterogeneity is increased. As can be seen in Fig. 1A, for $\sigma_h^2 = 0$ the potential $V(\bar{u})$ displays multiple optima, corresponding to the fixed points of Eq. (3). Further increases of $\sigma_h^2$ alters the shape of the potential $V(\bar{u})$, converging towards a purely quadratic form in which there is only one equilibrium. Hence, large heterogeneity yields trivialization.

C. Non-gradient systems

Conversely, general systems subscribing to Eq. (1) do
neither possess a well-defined mean field nor are they
gradient, being instead governed by random potentials²³.
We here summarize the theoretical framework proposed
in²³,²⁸, highlighting the relationship between the May-
Wigner transition threshold and the number of equilibria
of large-scale nonlinear systems such as those in Eq. (1)
in absence of a well-defined potential. For generality, we
focus on the well-studied case $\mu_w = 0$; we note however
that this condition may be relaxed and yield identical
results under some conditions¹⁴,²⁰. After the review of
major results on May-Wigner transitions in the following
paragraphs, sub-section II C 2, (Diversity Induced Trivial-
ization), provides new results on trivialization in random
networks with excitability heterogeneity.

1. Theory of May-Wigner transitions

One may indeed estimate the mean number of equilibria of Eq. (1), denoted by $\langle \mathcal{N} \rangle$, through the celebrated Kac-Rice relation$^{23-25,84}$

$$
\langle \mathcal{N} \rangle = \int_{\mathbb{R}^N} du \, \langle |\det(-I + J(\mathbf{u}))| \times \delta(L[\mathbf{u}] + Wf[\mathbf{u}]) \rangle \quad (6)
$$

where $\delta$ stands for the Dirac delta function, $\mathbf{J}$ the Jacobian associated with Eq. (1), $\mathbf{W}$ is the weight matrix with $(\mathbf{W})_{ij} = w_{ij}$, $(\mathbf{f}) = f_i$, $I$ is the $N$-dimensional identity matrix and $\langle \cdot \rangle$ corresponds to the average over realizations of the matrix $\mathbf{J}$. In the case where the Jacobian $\mathbf{J}$ is independent of $\mathbf{u}$, Eq. (6) simplifies further to$^{23}$

$$
\langle \mathcal{N} \rangle = \langle |\det(-\mathbf{I} + \mathbf{J})| \rangle \quad (7)
$$

Following the rationale detailed in²⁸, we know that for
spectral radii Γ < 1, Eq. (1) possesses a unique equilib-
rium denoted by **u**⁰. One further notices that

$$
\log |\det(-\mathbf{I} + \mathbf{J})| = \sum_{\lambda \in \Lambda} \log |\lambda - 1| \quad (8)
$$

where $\lambda \in \mathbb{C}$ are discrete eigenvalues and $\Lambda$ corresponds to the spectrum associated with the Jacobian $\mathbf{J}$. Dividing both sides of Eq. (8) by $N$ and taking the limit $N \to \infty$, one obtains

$$
\frac{1}{N} \sum_{\lambda \in \Lambda} \log |\lambda - 1| = \int_{\mathbb{C}} \log |s - 1| \rho(s) ds \quad (9)
$$

where $\rho(s)$ is the probability density of eigenvalues in the complex plane. For the system defined in Eq. (1), according to Girko's circular law of random matrices$^{76-78}$, the spectrum $\Lambda$ is uniformly distributed in the complex plane over a disk of radius $\Gamma$. For $s = re^{i\theta}$, $\rho(s) = \rho(r) = (\pi\Gamma^2)^{-1}$ for $r \le \Gamma$ and zero otherwise. Hence, one may write Eq. (9) as

$$
\int_{\mathbb{C}} \log |s-1| \rho(s) ds = \frac{1}{\pi \Gamma^2} \int_0^{2\pi} \int_0^{\Gamma} r \log |re^{i\theta} - 1| dr d\theta \quad (10)
$$

Noting that by residue calculus, one has

$$
\int_0^{2\pi} \log |re^{i\theta} - 1| d\theta = \begin{cases} 0 & 0 < r < 1 \\ 2\pi \log r & r \ge 1 \end{cases}
$$

And hence the right hand side of Eq. (10) simplifies to

$$
2\Gamma^{-2} \int_{1}^{\Gamma} r \log r \, dr = \log \Gamma + \frac{1}{2} \left( \frac{1}{\Gamma^2} - 1 \right)
$$

The expected number of equilibria $\langle \mathcal{N} \rangle$ is hence related,
to the first order, to the spectral radius $\Gamma$ by the following
well known$^{28}$ relationship

$$
\langle \mathcal{N} \rangle = \begin{cases} 1 & 0 < \Gamma < 1 \\ \exp[N \left( \log \Gamma + \frac{1}{2} \left( \frac{1}{\Gamma^2} - 1 \right) \right)] & \Gamma \ge 1 \end{cases} \quad (11)
$$

which confirms that whenever the spectral radius $\Gamma > 1$, an exponential increase in the number of equilibria accompanies the May-Wigner transition.

2. Diversity-Induced Trivialization

The quenching effect of excitability heterogeneity on
the number of equilibria for arbitrary non-gradient sys-
tems can be quantified through the same approach de-
scribed above. As per Eq. (11), $\langle \mathcal{N} \rangle$ depends, to the first
order, exclusively on the spectral radius $\Gamma$, which can be
computed explicitly in the presence of heterogeneity$^{65}$.
This allows us to revisit the above analysis and quantify
defined by

$$V(\bar{u}) = \frac{\bar{u}^2}{2} - \frac{\chi^o}{2\sqrt{\pi}\beta} \left( \gamma e^{\frac{-\beta^2(\mu_h - \bar{u})^2}{\gamma}} \right. \\ \left. + \sqrt{\pi}\beta \left( (\mu_h - \bar{u})\text{erf}\left[\frac{\beta(\mu_h - \bar{u})}{\sqrt{1 + 2\beta^2\sigma_h^2}}\right] + \bar{u} \right) \right) \quad (5)$$

Figure 1 depicts trivialization that occurs in Eq. (3) once excitability heterogeneity is increased. As can be seen in Fig. 1A, for $\sigma_h^2 = 0$ the potential $V(\bar{u})$ displays multiple optima, corresponding to the fixed points of Eq. (3). Further increases of $\sigma_h^2$ alters the shape of the potential $V(\bar{u})$, converging towards a purely quadratic form in which there is only one equilibrium. Hence, large heterogeneity yields trivialization.

### C. Non-gradient systems

Conversely, general systems subscribing to Eq. (1) do neither possess a well-defined mean field nor are they gradient, being instead governed by random potentials$^{23}$. We here summarize the theoretical framework proposed in$^{23,28}$, highlighting the relationship between the May-Wigner transition threshold and the number of equilibria of large-scale nonlinear systems such as those in Eq. (1) in absence of a well-defined potential. For generality, we focus on the well-studied case $\mu_w = 0$; we note however that this condition may be relaxed and yield identical results under some conditions$^{14,20}$. After the review of major results on May-Wigner transitions in the following paragraphs, sub-section II C 2, (Diversity Induced Trivialization), provides new results on trivialization in random networks with excitability heterogeneity.

#### 1. Theory of May-Wigner transitions

One may indeed estimate the mean number of equilibria of Eq. (1), denoted by $\langle \mathcal{N} \rangle$, through the celebrated Kac-Rice relation$^{23-25,84}$

$$\langle \mathcal{N} \rangle = \int_{\mathbb{R}^N} d\mathbf{u} \, \langle |\det(-\mathbb{I} + J(\mathbf{u}))| \times \delta(L[\mathbf{u}] + Wf[\mathbf{u}]) \rangle \quad (6)$$

where $\delta$ stands for the Dirac delta function, $J$ the Jacobian associated with Eq. (1), $W$ is the weight matrix with $(W)_{ij} = w_{ij}$, $(f) = f_i$, $\mathbb{I}$ is the $N$-dimensional identity matrix and $\langle \cdot \rangle$ corresponds to the average over realizations of the matrix $J$. In the case where the Jacobian $J$ is independent of $\mathbf{u}$, Eq. (6) simplifies further to$^{23}$

$$\langle \mathcal{N} \rangle = \langle |\det(-\mathbb{I} + J)| \rangle \quad (7)$$

Following the rationale detailed in$^{28}$, we know that for spectral radii $\Gamma < 1$, Eq. (1) possesses a unique equilibrium denoted by $\mathbf{u}^o$. One further notices that

$$\log |\det(-\mathbb{I} + J)| = \sum_{\lambda \in \Lambda} \log |\lambda - 1| \quad (8)$$

where $\lambda \in \mathbb{C}$ are discrete eigenvalues and $\Lambda$ corresponds to the spectrum associated with the Jacobian $J$. Dividing both sides of Eq. (8) by $N$ and taking the limit $N \to \infty$, one obtains

$$\frac{1}{N} \sum_{\lambda \in \Lambda} \log |\lambda - 1| = \int_{\mathbb{C}} \log |s - 1| \rho(s) ds \quad (9)$$

where $\rho(s)$ is the probability density of eigenvalues in the complex plane. For the system defined in Eq. (1), according to Girko's circular law of random matrices$^{76-78}$, the spectrum $\Lambda$ is uniformly distributed in the complex plane over a disk of radius $\Gamma$. For $s = re^{i\theta}$, $\rho(s) = \rho(r) = (\pi\Gamma^2)^{-1}$ for $r \le \Gamma$ and zero otherwise. Hence, one may write Eq. (9) as

$$\int_{\mathbb{C}} \log |s - 1| \rho(s) ds = \frac{1}{\pi \Gamma^2} \int_0^{2\pi} \int_0^{\Gamma} r \log |re^{i\theta} - 1| dr d\theta \quad (10)$$

Noting that by residue calculus, one has

$$\int_{0}^{2\pi} \log |re^{i\theta} - 1| d\theta = \begin{cases} 0 & 0 < r < 1 \\ 2\pi \log r & r \ge 1 \end{cases}$$

And hence the right hand side of Eq. (10) simplifies to

$$2\Gamma^{-2} \int_{1}^{\Gamma} r \log r dr = \log \Gamma + \frac{1}{2} \left( \frac{1}{\Gamma^2} - 1 \right)$$

The expected number of equilibria $\langle \mathcal{N} \rangle$ is hence related, to the first order, to the spectral radius $\Gamma$ by the following well known$^{28}$ relationship

$$\langle \mathcal{N} \rangle = \begin{cases} 1 & 0 < \Gamma < 1 \\ \exp[N \left( \log \Gamma + \frac{1}{2} \left( \frac{1}{\Gamma^2} - 1 \right) \right)] & \Gamma \ge 1 \end{cases} \quad (11)$$

which confirms that whenever the spectral radius $\Gamma > 1$, an exponential increase in the number of equilibria accompanies the May-Wigner transition.

#### 2. Diversity-Induced Trivialization

The quenching effect of excitability heterogeneity on the number of equilibria for arbitrary non-gradient systems can be quantified through the same approach described above. As per Eq. (11), $\langle \mathcal{N} \rangle$ depends, to the first order, exclusively on the spectral radius $\Gamma$, which can be computed explicitly in the presence of heterogeneity$^{65}$. This allows us to revisit the above analysis and quantify
how $\langle \mathcal{N} \rangle$ scales with the dispersion of intrinsic neuronal properties.

One may note that below the May-Wigner transition ($\Gamma < 1$), as per the Banach fixed point theorem, a unique stable fixed point $\mathbf{u}^o = [u_1^o, u_2^o, ..., u_N^o]$ exists, satisfying the set of implicit equations

$$u_i^o = \sum_{j=1}^{N} w_{ij} f_j[u_j^o] \quad (12)$$

Now let us consider neuron-specific excitability thresholds $h_j \neq 0$ in Eq. (1) where $h_j$ obey the Gaussian probability density function $\rho(h)$ with vanishing mean and non-vanishing variance $\sigma_h^2$. Then the resulting unique fixed point $\mathbf{u}^o$ is not trivial and distributed with a probability density function $\rho(\mathbf{u}^o)$. This function's first two moments can be easily computed, given a few simplifying assumptions. Indeed, assuming statistical independence between firing rates, synaptic weights and heterogeneity, and whenever $N$ is sufficiently large, the mean over the nodes in Eq. (12) may now be approximated by the following statistical average

$$\mu_{u^o} \equiv \langle u_i^o \rangle_N = \left\langle \sum_{j=1}^{N} w_{ij} f_j[u_j^o] \right\rangle_N \\ = \bar{\mu}_w \int_{-\infty}^{\infty} f[u^o] \rho(u^o) du^o = 0 \quad (13)$$

with $\bar{\mu}_w = N\mu_w$ and by independence of connectivity weights and fixed point coordinates. Equation (13) is an implicit equation for $\mu_{u^o}$ and holds true since $f[u^o]$ and $\rho(u^o)$ are anti-symmetric and symmetric to the origin, respectively.

The variance $\sigma_{\bar{u}^o}^2$ can be computed similarly$^{65}$

$$\sigma_{\bar{u}^o}^2 = \langle (u_i^o)^2 \rangle_N = \left\langle \sum_{j,k=1}^{N} w_{ij} w_{ik} f_j[u_j^o] f_k[u_k^o] \right\rangle_N \\ = \frac{N \sigma_W^2}{4} \left( 1 - \frac{2}{\sqrt{4 + \beta^2 \pi^2 \sigma_h^2}} \right) . \quad (14)$$

We proceed by first observing that elements of the Jacobian $J$ associated with Eq. (12) can be written as

$$\begin{aligned} J_{ij} &= -\delta_{ij} + w_{ij} \frac{\partial f_j}{\partial u_j} |_{u_j^o} \\ &= -\delta_{ij} + \frac{w_{ij} \beta}{\sqrt{\pi}} e^{-\beta^2(u_j^o+h_j)^2} \end{aligned} \quad (15)$$

Substituting two normally distributed and independent random variables $s$ and $q$ whose statistics are identical to those of $u_i^o$ and $h_j$ i.e. $\rho(s) \to \rho(\mathbf{u}^o)$ and $\rho(q) \to \rho(h)$, one obtains a surrogate yet statistically identical Jacobian

$\tilde{J}$ i.e.

$$\begin{aligned} \tilde{J}_{ij} &= J_{ij} + \delta_{ij} \\ &= \tilde{w}_{ij} \end{aligned} \quad (16)$$

where we have used the change of variable, effectively re-scaling the synaptic weights,

$$\tilde{w}_{ij} = \frac{w_{ij} \beta}{\sqrt{\pi}} e^{-\beta^2 (s+q)^2}$$

Through this substitution, the calculation of the spectral radius follows$^{9,10,14,15}$

$$\Gamma = \sqrt{(N-1)\rho\sigma_{\bar{w}}^2}, \quad (17)$$

i.e. the fixed point $\mathbf{u}^o$ is stable if $\Gamma < 1$.

In addition, whenever $N$ is sufficiently large, the variance of the re-scaled synaptic weights becomes,

$$\begin{aligned} \sigma_{\bar{W}}^2 &= \sigma_{\bar{W}}^2 \mathrm{Var} \left[ \frac{\beta}{\sqrt{\pi}} e^{-\beta^2 (s+q)^2} \right] \\ &= \frac{\sigma_{\bar{W}}^2 \beta^2}{\pi} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} e^{-\beta^2 (s+q)^2} \rho(s) \rho(q) ds dq \\ &= \frac{\sigma_{\bar{W}}^2 \beta^2}{\pi \sqrt{\gamma}} \end{aligned} \quad (18)$$

with

$$\gamma = 1 + 4\beta^2 (\sigma_{\bar{u}^o}^2 + \sigma_h^2) .$$

Hence the spectral radius reads

$$\Gamma = \sqrt{\frac{\sigma_{\bar{W}}^2 (N-1) \rho \beta^2}{\pi \sqrt{\gamma}}} . \quad (19)$$

utilizing Eq. (14).

We emphasize that Eq. (19) applies for Erdös-Renyi connectivity matrices, however, our results can easily be extended to other topologies, such as multi-modal, scale-free or cascade models, through a proper re-scaling of the spectral radius$^{22,27}$.

By virtue of the statistical equivalence of the Jacobian $J$ and $\tilde{J}$ and their independence with respect to $\mathbf{u}^o$, the expected number of fixed points may hence be approximated through Eq. (7)

$$\langle \mathcal{N} \rangle (\sigma_h^2) \approx \langle |\det(-I + \tilde{J})| \rangle$$

and Eq. (11) follows directly. We may hence substitute Eq. (19) in Eq. (11) to obtain an explicit expression for the expected number of equilibria $\langle \mathcal{N} \rangle$ as a function of the variance in excitability thresholds. Specifically, we may expand Eq. (11) at the stability threshold $\Gamma \sim 1$ and obtain

$$\langle \mathcal{N} \rangle (\sigma_h^2) \approx 1 + N \left( \sqrt{\frac{\sigma_{\bar{W}}^2 (N-1) \rho \beta^2}{\pi \sqrt{\gamma}}} - 1 \right)^2 . \quad (20)$$
how $\langle \mathcal{N} \rangle$ scales with the dispersion of intrinsic neuronal properties.

One may note that below the May-Wigner transition ($\Gamma < 1$), as per the Banach fixed point theorem, a unique stable fixed point $\mathbf{u}^o = [u_1^o, u_2^o, ..., u_N^o]$ exists, satisfying the set of implicit equations

$$u_i^o = \sum_{j=1}^{N} w_{ij} f_j[u_j^o] \quad (12)$$

Now let us consider neuron-specific excitability thresholds $h_j \neq 0$ in Eq. (1) where $h_j$ obey the Gaussian probability density function $\rho(h)$ with vanishing mean and non-vanishing variance $\sigma_h^2$. Then the resulting unique fixed point $\mathbf{u}^o$ is not trivial and distributed with a probability density function $\rho(\mathbf{u}^o)$. This function's first two moments can be easily computed, given a few simplifying assumptions. Indeed, assuming statistical independence between firing rates, synaptic weights and heterogeneity, and whenever $N$ is sufficiently large, the mean over the nodes in Eq. (12) may now be approximated by the following statistical average

$$\mu_{u^o} \equiv \langle u_i^o \rangle_N = \left\langle \sum_{j=1}^{N} w_{ij} f_j[u_j^o] \right\rangle_N \\ = \bar{\mu}_w \int_{-\infty}^{\infty} f[u^o] \rho(u^o) du^o = 0 \quad (13)$$

with $\bar{\mu}_w = N\mu_w$ and by independence of connectivity weights and fixed point coordinates. Equation (13) is an implicit equation for $\mu_{u^o}$ and holds true since $f[u^o]$ and $\rho(u^o)$ are anti-symmetric and symmetric to the origin, respectively.

The variance $\sigma_{\bar{u}^o}^2$ can be computed similarly$^{65}$

$$\sigma_{\bar{u}^o}^2 = \langle (u_i^o)^2 \rangle_N = \left\langle \sum_{j,k=1}^{N} w_{ij} w_{ik} f_j[u_j^o] f_k[u_k^o] \right\rangle_N \\ = \frac{N \sigma_W^2}{4} \left( 1 - \frac{2}{\sqrt{4 + \beta^2 \pi^2 \sigma_h^2}} \right) . \quad (14)$$

We proceed by first observing that elements of the Jacobian $J$ associated with Eq. (12) can be written as

$$\begin{aligned} J_{ij} &= -\delta_{ij} + w_{ij} \frac{\partial f_j}{\partial u_j} |_{u_j^o} \\ &= -\delta_{ij} + \frac{w_{ij} \beta}{\sqrt{\pi}} e^{-\beta^2(u_j^o+h_j)^2} \end{aligned} \quad (15)$$

Substituting two normally distributed and independent random variables $s$ and $q$ whose statistics are identical to those of $u_i^o$ and $h_j$ i.e. $\rho(s) \to \rho(\mathbf{u}^o)$ and $\rho(q) \to \rho(h)$, one obtains a surrogate yet statistically identical Jacobian

$\tilde{J}$ i.e.

$$\begin{aligned} \tilde{J}_{ij} &= J_{ij} + \delta_{ij} \\ &= \tilde{w}_{ij} \end{aligned} \quad (16)$$

where we have used the change of variable, effectively re-scaling the synaptic weights,

$$\tilde{w}_{ij} = \frac{w_{ij} \beta}{\sqrt{\pi}} e^{-\beta^2 (s+q)^2}$$

Through this substitution, the calculation of the spectral radius follows$^{9,10,14,15}$

$$\Gamma = \sqrt{(N-1)\rho\sigma_{\bar{w}}^2}, \quad (17)$$

i.e. the fixed point $\mathbf{u}^o$ is stable if $\Gamma < 1$.

In addition, whenever $N$ is sufficiently large, the variance of the re-scaled synaptic weights becomes,

$$\begin{aligned} \sigma_{\bar{w}}^2 &= \sigma_W^2 \text{Var} \left[ \frac{\beta}{\sqrt{\pi}} e^{-\beta^2 (s+q)^2} \right] \\ &= \frac{\sigma_W^2 \beta^2}{\pi} \int_{-\infty}^{\infty} \int_{-\infty}^{\infty} e^{-\beta^2 (s+q)^2} \rho(s) \rho(q) ds dq \\ &= \frac{\sigma_W^2 \beta^2}{\pi \sqrt{\gamma}} \end{aligned} \quad (18)$$

with

$$\gamma = 1 + 4\beta^2 (\sigma_{\bar{u}^o}^2 + \sigma_h^2) .$$

Hence the spectral radius reads

$$\Gamma = \sqrt{\frac{\sigma_W^2 (N-1) \rho \beta^2}{\pi \sqrt{\gamma}}} . \quad (19)$$

utilizing Eq. (14).

We emphasize that Eq. (19) applies for Erdös-Renyi connectivity matrices, however, our results can easily be extended to other topologies, such as multi-modal, scale-free or cascade models, through a proper re-scaling of the spectral radius$^{22,27}$.

By virtue of the statistical equivalence of the Jacobian $J$ and $\tilde{J}$ and their independence with respect to $\mathbf{u}^o$, the expected number of fixed points may hence be approximated through Eq. (7)

$$\langle \mathcal{N} \rangle (\sigma_h^2) \approx \langle |\det(-I + \tilde{J})| \rangle$$

and Eq. (11) follows directly. We may hence substitute Eq. (19) in Eq. (11) to obtain an explicit expression for the expected number of equilibria $\langle \mathcal{N} \rangle$ as a function of the variance in excitability thresholds. Specifically, we may expand Eq. (11) at the stability threshold $\Gamma \sim 1$ and obtain

$$\langle \mathcal{N} \rangle (\sigma_h^2) \approx 1 + N \left( \sqrt{\frac{\sigma_W^2 (N-1) \rho \beta^2}{\pi \sqrt{\gamma}}} - 1 \right)^2 . \quad (20)$$
Obviously, the mean number of fixed points $\langle N \rangle$ is a decaying function with respect to the excitability heterogeneity variance $\sigma_h^2$. This confirms that trivialization occurs whenever intrinsic excitability increases its variance. This translates into an effective shift in the May-Wigner transition threshold occurring at $\Gamma = 1$ for a critical heterogeneity level$^{54,56,62}$ $\sigma_{h,c}^2$. For $\sigma_h^2 \ll 1$, one may expand Eq. (14) to first order in $\sigma_h^2$ to obtain an explicit expression for the critical heterogeneity for trivialization, i.e.

$$ \sigma_{h,c}^2 \approx \frac{8}{\pi^2 \beta^2} \left( \frac{N^2 \beta^4 \sigma_W^4 - \pi^2}{N \sigma_W^2 \pi^2 \beta^2 + 32} \right) \quad (21) $$

It has been shown$^{65}$ that dispersion of excitability heterogeneity (i.e. $\sigma_h^2 > 0$) effectively contracts the spectral radius, stabilizing the dynamics and decreasing the probability of state transitions. This phenomenon translates here directly to trivialization: as can be seen in Fig. 2, both numerical and analytical calculations jointly predict a decrease in the number of expected equilibria $\langle N \rangle$ as the dispersion of excitability thresholds increases. Such diversity seems to limit the exponential explosion of fixed points and appears to delay the May-Wigner transition generically observed whenever $\Gamma \ge 1$.

### III. APPLICATIONS

#### A. Resilience

We have demonstrated in the previous section, that excitability heterogeneity induces trivialization. Since trivialization is strongly related to the resilience of a system to perturbations, the current section aims to demonstrate enhanced resilience by stronger excitability heterogeneity. We assume that a system is in its stationary state and receives a short strong perturbation impulse $I(t)$, that repels the system from its stationary state. After the perturbation impulse, the system evolves according to its unperturbed dynamics and relaxes to a new stationary state. We state that the system is resilient if the new stationary state is the same as the original state, where the system was located before the perturbation. We quantify the degree of resilience, and compute the distance $D$ of the newly approached state to the original state. If $D \approx 0$, then the system is resilient, whereas $D > 0$ indicates lower resilience.

In more detail, we assume model (1) with $L[u] = -u$ and $\mu_w = 1.0$, $\sigma_w^2 = 0.05$, $\mu_h = 1.0$ and $\sigma_h^2 > 0$. The distance reflecting the degree of resilience is defined as

$$ D = \frac{1}{N} \sum_{i=1}^N |u_i^{\text{original}} - u_i^{\text{perturbed}}| . \quad (22) $$

We simulated Eq. (1) with initial condition $u_i(t = 0) = \xi_i$, $i=1, \dots, N$, where $\xi_i \in [0; 1]$ has been uniform distributed. Temporal integration with time step $\Delta = 0.05$

yielded a stationary state if $\sum_{i=1}^N |\Delta u_i|/N < 10^{-3}$ with $\Delta u_i = -u_i + \sum_{j=1}^N w_{ij} f_j[u_j]$. Then a random uniform-distributed instantaneous perturbation $Y_i \in [-500; 500]$ is applied at each network node $i$, after which the system approaches the new stationary state $u_i^{\text{perturbed}}$. This new state is reached if $\sum_{i=1}^N |\Delta u_i|/N < 10^{-4}$. For each excitability variance $\sigma_h^2$, we have run $T = 400$ simulation trials, where each trial draws new random connections $w_{ij}$ and new excitability values $h_i$ from their corresponding distributions. This yields a value of $D$ for each $\sigma_h^2$ and we have computed this value $T$ times. Figure 3 shows all obtained values of $D$ for difference values of the excitability heterogeneity variance $\sigma_h^2$. We observe that systems with weak heterogeneity (low values of $\sigma_h^2$) exhibit different non-zero values of $D$ and thus a low degree of resilience, whereas increasing heterogeneity yields better and better resilience.

#### B. Neural development and learning

The current work is strongly motivated by the role of intrinsic diversity during development, through which synaptic connectivity motifs become increasingly dense and complex. What maintains neural circuit stability? Indeed, during development and learning, the number of synapses and their effective weights change, notably through long-term synaptic potentiation (LTP) and depression (LTD)$^{85}$. Such increase in complexity, parameterized by gradual increase in connection probability and/or density of synaptic connections and/or changes to the number of neurons and synapses, was repeatedly found to foster instability by increasing the risk of May-Wigner type transitions$^{8–12,14,18,19,21,22,86}$. Similarly, synaptic plasticity models in recurrent neural networks commonly face unstable regimes which need to be compensated by/stabilized through various feedback processes$^{33,87–89}$. While many neurophysiological processes, such as synaptic inhibition, adaptation, modularity, and circuit motifs are certainly involved in stabilizing neural circuits’ function, we believe that intrinsic excitability heterogeneity is also involved through trivialization. To emulate such a developmental/learning process for networks whose dynamics obey Eq. (1), we considered a population of $N$ initially isolated neurons, whose mutual connectivity would increase gradually in time (cf. Fig. 4). To implement this, we simulated the solutions of Eq. (1) while the connection probability $\rho$ linearly increases from 0 to 5% over an extended time scale (Fig. 4A). Note that while synaptic weights would also be expected to increase during development/learning, the influence of this parameter on stability would be identical (cf. Eq. (19)). Hence we focused on the connection probability for simplicity. Results presented in Fig. 4B confirm the destabilizing effect of increase in connectivity and the resurgence of a May-Wigner instability: the network dynamics reflects a
Obviously, the mean number of fixed points $\langle N \rangle$ is a decaying function with respect to the excitability heterogeneity variance $\sigma_h^2$. This confirms that trivialization occurs whenever intrinsic excitability increases its variance. This translates into an effective shift in the May-Wigner transition threshold occurring at $\Gamma = 1$ for a critical heterogeneity level$^{54,56,62}$ $\sigma_{h,c}^2$. For $\sigma_h^2 \ll 1$, one may expand Eq. (14) to first order in $\sigma_h^2$ to obtain an explicit expression for the critical heterogeneity for trivialization, i.e.

$$ \sigma_{h,c}^2 \approx \frac{8}{\pi^2 \beta^2} \left( \frac{N^2 \beta^4 \sigma_W^4 - \pi^2}{N \sigma_W^2 \pi^2 \beta^2 + 32} \right) \quad (21) $$

It has been shown$^{65}$ that dispersion of excitability heterogeneity (i.e. $\sigma_h^2 > 0$) effectively contracts the spectral radius, stabilizing the dynamics and decreasing the probability of state transitions. This phenomenon translates here directly to trivialization: as can be seen in Fig. 2, both numerical and analytical calculations jointly predict a decrease in the number of expected equilibria $\langle N \rangle$ as the dispersion of excitability thresholds increases. Such diversity seems to limit the exponential explosion of fixed points and appears to delay the May-Wigner transition generically observed whenever $\Gamma \ge 1$.

### III. APPLICATIONS

#### A. Resilience

We have demonstrated in the previous section, that excitability heterogeneity induces trivialization. Since trivialization is strongly related to the resilience of a system to perturbations, the current section aims to demonstrate enhanced resilience by stronger excitability heterogeneity. We assume that a system is in its stationary state and receives a short strong perturbation impulse $I(t)$, that repels the system from its stationary state. After the perturbation impulse, the system evolves according to its unperturbed dynamics and relaxes to a new stationary state. We state that the system is resilient if the new stationary state is the same as the original state, where the system was located before the perturbation. We quantify the degree of resilience, and compute the distance $D$ of the newly approached state to the original state. If $D \approx 0$, then the system is resilient, whereas $D > 0$ indicates lower resilience.

In more detail, we assume model (1) with $L[u] = -u$ and $\mu_w = 1.0$, $\sigma_w^2 = 0.05$, $\mu_h = 1.0$ and $\sigma_h^2 > 0$. The distance reflecting the degree of resilience is defined as

$$ D = \frac{1}{N} \sum_{i=1}^N |u_i^{\text{original}} - u_i^{\text{perturbed}}| . \quad (22) $$

We simulated Eq. (1) with initial condition $u_i(t = 0) = \xi_i$, $i=1, \dots, N$, where $\xi_i \in [0; 1]$ has been uniform distributed. Temporal integration with time step $\Delta = 0.05$

yielded a stationary state if $\sum_{i=1}^N |\Delta u_i|/N < 10^{-3}$ with $\Delta u_i = -u_i + \sum_{j=1}^N w_{ij} f_j[u_j]$. Then a random uniform-distributed instantaneous perturbation $Y_i \in [-500; 500]$ is applied at each network node $i$, after which the system approaches the new stationary state $u_i^{\text{perturbed}}$. This new state is reached if $\sum_{i=1}^N |\Delta u_i|/N < 10^{-4}$. For each excitability variance $\sigma_h^2$, we have run $T = 400$ simulation trials, where each trial draws new random connections $w_{ij}$ and new excitability values $h_i$ from their corresponding distributions. This yields a value of $D$ for each $\sigma_h^2$ and we have computed this value $T$ times. Figure 3 shows all obtained values of $D$ for difference values of the excitability heterogeneity variance $\sigma_h^2$. We observe that systems with weak heterogeneity (low values of $\sigma_h^2$) exhibit different non-zero values of $D$ and thus a low degree of resilience, whereas increasing heterogeneity yields better and better resilience.

#### B. Neural development and learning

The current work is strongly motivated by the role of intrinsic diversity during development, through which synaptic connectivity motifs become increasingly dense and complex. What maintains neural circuit stability? Indeed, during development and learning, the number of synapses and their effective weights change, notably through long-term synaptic potentiation (LTP) and depression (LTD)$^{85}$. Such increase in complexity, parameterized by gradual increase in connection probability and/or density of synaptic connections and/or changes to the number of neurons and synapses, was repeatedly found to foster instability by increasing the risk of May-Wigner type transitions$^{8–12,14,18,19,21,22,86}$. Similarly, synaptic plasticity models in recurrent neural networks commonly face unstable regimes which need to be compensated by/stabilized through various feedback processes$^{33,87–89}$. While many neurophysiological processes, such as synaptic inhibition, adaptation, modularity, and circuit motifs are certainly involved in stabilizing neural circuits’ function, we believe that intrinsic excitability heterogeneity is also involved through trivialization. To emulate such a developmental/learning process for networks whose dynamics obey Eq. (1), we considered a population of $N$ initially isolated neurons, whose mutual connectivity would increase gradually in time (cf. Fig. 4). To implement this, we simulated the solutions of Eq. (1) while the connection probability $\rho$ linearly increases from 0 to 5% over an extended time scale (Fig. 4A). Note that while synaptic weights would also be expected to increase during development/learning, the influence of this parameter on stability would be identical (cf. Eq. (19)). Hence we focused on the connection probability for simplicity. Results presented in Fig. 4B confirm the destabilizing effect of increase in connectivity and the resurgence of a May-Wigner instability: the network dynamics reflects a
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Example of diversity-induced trivialization in a random neural network A. Response of a homogeneous network, in which all neurons have the same excitability profile (i.e. $h_j = 0 \forall j$). Chaotic dynamics can be observed, indicative that the system resides beyond the May-Wigner transition threshold B. Associated eigenvalues in the complex plane (black dots) alongside the theoretically derived spectral radius (19) plotted around $\lambda = -1$ (blue circle). Some eigenvalues occupy the right hand side of the imaginary axis, confirming instability. C-D. Dynamics and eigenspectrum of the same network whenever heterogeneity is increased to $\sigma_h^2 = 0.0001$ and; E-F $\sigma_h^2 = 0.001$. Heterogeneity stabilizes the system by contracting the spectral radius in the complex plane, yielding stable dynamics. Blue dots in C,E illustrate different values of $h_j$. G Spectral radius $\Gamma$ as a function of excitability heterogeneity. As $\sigma_h^2$ increases and passes the critical point $\sigma_{h,c}^2$ (vertical dashed line), the May-Wigner transition threshold (i.e. $\Gamma = 1$; horizontal line) is crossed, stabilizing the dynamics (cf Panels A-F). Numerical results averaged over 50 independent realizations of the network connectivity (black dotted line; grey shaded area represents $\pm$ SD) are plotted alongside theoretical predictions (solid blue line). H Expected number of equilibria $\langle N \rangle$ as a function of $\sigma_h^2$, revealing trivialization and a shift in the May-Wigner instability threshold induced by excitability heterogeneity. Numerical results (black dotted line) and theoretical predictions (solid blue line) coincide at and beyond the critical point, but diverge otherwise due to error propagation through an exponential (see Eq. (11)). I Variance of the steady state probability density function as a function of excitability heterogeneity, confirming that $\sigma_h^2$ does increase the dispersion of the steady state, while showing good alignment between numerical (black dotted line; grey shaded area represents $\pm$ SD) and theoretical (solid blue line) calculations. J Variance of the dynamics for the individual neurons (i.e. $u_j$; grey shaded lines) across various levels of excitability heterogeneity, shown alongside their average over independent trials (black dotted line). As the system transits from chaotic to stable dynamics by crossing the critical threshold $\sigma_{h,c}^2$ (cf. panels A-F), the system activity variance decreases. In other words, increasing the excitability heterogeneity reduces the system fluctuations and increases system order. Parameters are $N = 100$, $\rho = 0.05$, $g = 0.8$, $\mu_e = 0.005$, $\beta = 25$, $\mu_i = \frac{g}{g-1}\mu_e$, $\sigma_{W,e}^2 = \sigma_{W,i}^2 = 0.0015$. Whenever applicable, averages were computed over 50 independent trials of duration $T = 200$ [a.u.].</figcaption></figure>

highly unstable, multi-modal oscillatory state. Introducing excitability heterogeneity in Fig. 4C was in contrast found to stabilize the dynamics by dispersing the neurons' activity. Spectral radii calculations affirmed these numerical observations (Fig. 4D) and further revealed

that heterogeneity delays the May-Wigner transition and its associated exponential increase in equilibria (Fig. 4E). Trivialization was further found to scale with the level of excitability heterogeneity (i.e. $\sigma_h^2$): as can be seen in Fig. 4F-G, increased level of heterogeneity decreases
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Example of diversity-induced trivialization in a random neural network A. Response of a homogeneous network, in which all neurons have the same excitability profile (i.e. $h_j = 0 \forall j$). Chaotic dynamics can be observed, indicative that the system resides beyond the May-Wigner transition threshold B. Associated eigenvalues in the complex plane (black dots) alongside the theoretically derived spectral radius (19) plotted around $\lambda = -1$ (blue circle). Some eigenvalues occupy the right hand side of the imaginary axis, confirming instability. C-D. Dynamics and eigenspectrum of the same network whenever heterogeneity is increased to $\sigma_h^2 = 0.0001$ and; E-F $\sigma_h^2 = 0.001$. Heterogeneity stabilizes the system by contracting the spectral radius in the complex plane, yielding stable dynamics. Blue dots in C,E illustrate different values of $h_j$. G Spectral radius $\Gamma$ as a function of excitability heterogeneity. As $\sigma_h^2$ increases and passes the critical point $\sigma_{h,c}^2$ (vertical dashed line), the May-Wigner transition threshold (i.e. $\Gamma = 1$; horizontal line) is crossed, stabilizing the dynamics (cf Panels A-F). Numerical results averaged over 50 independent realizations of the network connectivity (black dotted line; grey shaded area represents $\pm$ SD) are plotted alongside theoretical predictions (solid blue line). H Expected number of equilibria $\langle N \rangle$ as a function of $\sigma_h^2$, revealing trivialization and a shift in the May-Wigner instability threshold induced by excitability heterogeneity. Numerical results (black dotted line) and theoretical predictions (solid blue line) coincide at and beyond the critical point, but diverge otherwise due to error propagation through an exponential (see Eq. (11)). I Variance of the steady state probability density function as a function of excitability heterogeneity, confirming that $\sigma_h^2$ does increase the dispersion of the steady state, while showing good alignment between numerical (black dotted line; grey shaded area represents $\pm$ SD) and theoretical (solid blue line) calculations. J Variance of the dynamics for the individual neurons (i.e. $u_j$; grey shaded lines) across various levels of excitability heterogeneity, shown alongside their average over independent trials (black dotted line). As the system transits from chaotic to stable dynamics by crossing the critical threshold $\sigma_{h,c}^2$ (cf. panels A-F), the system activity variance decreases. In other words, increasing the excitability heterogeneity reduces the system fluctuations and increases system order. Parameters are $N = 100$, $\rho = 0.05$, $g = 0.8$, $\mu_e = 0.005$, $\beta = 25$, $\mu_i = \frac{g}{g-1}\mu_e$, $\sigma_{W,e}^2 = \sigma_{W,i}^2 = 0.0015$. Whenever applicable, averages were computed over 50 independent trials of duration $T = 200$ [a.u.].</figcaption></figure>

highly unstable, multi-modal oscillatory state. Introducing excitability heterogeneity in Fig. 4C was in contrast found to stabilize the dynamics by dispersing the neurons' activity. Spectral radii calculations affirmed these numerical observations (Fig. 4D) and further revealed

that heterogeneity delays the May-Wigner transition and its associated exponential increase in equilibria (Fig. 4E). Trivialization was further found to scale with the level of excitability heterogeneity (i.e. $\sigma_h^2$): as can be seen in Fig. 4F-G, increased level of heterogeneity decreases
<figure><img src="image_5.png" /><figcaption>**Figure 3. Resilience by excitability heterogeneity** The distance *D* taken from Eq. (22) reflects the degree of the system's resilience: the smaller *D*, the closer are the original and perturbed state and the more resilient is the system towards perturbations. We observe that strong heterogeneity renders the system more resilient to perturbations. For each value of $\sigma_h^2$, 400 values of *D* are plotted to correspond with the same number of simulation trials.</figcaption></figure>

the spectral radius and further limits the number of fixed points. Taken together, these findings suggest that heterogeneity quenches the dynamical repertoire of non-linear networks as their complexity increases to maintain their resilience, suggesting that neural diversity plays a role in the persistence of brain function.

## IV. DISCUSSION

Excitability heterogeneity affects stability of random networks and hence can induce trivialization and enhance resilience. This is shown theoretically in Eq. (5) in gradient systems and in the corresponding application shown in Fig. 3. In non-gradient systems, Equation (11) confirms that an exponential increase in fixed points occurs whenever the May-Wigner instability threshold is met i.e. $\Gamma \ge 1$$^{23-25,28}$. An important corollary of this equation is that the number of equilibria scales with the spectral radius, whose magnitude has been shown to increase with a variety of topological traits, such as increased dimensionality$^{9,18,19}$, strength of connections$^{10-12,14,18,21}$, motifs richness$^{10-12,14,18,21}$, connectance$^{8,9,14,18,21}$ and degree heterogeneity$^{22}$. Interestingly, our results as well as others$^{65}$ indicate that increased intrinsic diversity in cellular excitability represents an exception, contracting the spectral radius and hence promoting stability through a reduction of the number of equilibria. This suggests that variability amongst cells of a network - as opposed to heterogeneity in their connectivity - promotes resilience, by reducing the probability of multistability and bifurcations.

It is interesting to conjecture on the potential role of trivialization in neural circuits' dynamics and function. Intuitively, increasing the number of control parameters by enabling cells with different excitability profiles should enrich a system's phase space – however, our results suggest

the opposite. Numerous experimental$^{2,44-47,49-52}$ and theoretical studies$^{2,53-59}$ have examined the role of cell-to-cell heterogeneity in brain dynamics and function. For instance, several neural population models considered heterogeneous single cell properties$^{90-92}$. Collectively, these indicate that intrinsic diversity may both enhance resilience$^{47,48,52,63,64}$ and optimize coding performance$^{54,56,58,59,62}$. Our complementary results show that variations in cellular excitability represents a stabilizing mechanism, quenching the dynamical repertoire of balanced neural networks. This surprising feature arises notably through the effective linearization of input-output response functions$^{2,23,93}$, promoting graded response to perturbations and subsequent increasingly linear relaxation towards a limited number of stable states.

Multiple mechanisms have been proposed to regularize the dynamics of complex networks beyond heterogeneity. Noise and disorder - ubiquitous in neural systems$^{94}$ - represent important drivers of regularization, where random-like perturbations may paradoxically yield highly predictable dynamics over a wide range of physical$^{67-69}$ and biological$^{70,71}$ systems. This form of complexity collapse$^{95}$ has important implications for neural circuit function. These suggest a trade-off between topological complexity – which has been extensively studied – and intrinsic diversity, balancing coding efficiency and resilience. As such, trivialization resulting from either noise, disorder or heterogeneity, should be considered as a control mechanism, tuning stability towards an optimal dynamical state$^{96-100}$. For instance, destabilizing changes in connectivity during development and learning may be compensated for through adaptive increase in the level of fluctuations, biological diversity or both, to optimize neural circuit performance.

Intrinsic plasticity$^{33-37}$ refers to a non-synaptic form of plasticity through which individual cells' biophysical characteristics and/or excitability change in experience and stimuli-dependent ways, notably through adjustments in ionic conductances$^{47,52,63,101}$. While speculative, our results nonetheless suggest that intrinsic plasticity could play an important role in stabilizing neural dynamics and function, compensating for the notoriously destabilizing effect of synaptic plasticity$^{33,87-89}$.

## DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## ACKNOWLEDGMENTS

We thank the National Sciences and Engineering Research Council of Canada (NSERC Grants RGPIN-2017-
<figure><img src="image_4.png" /><figcaption>**Figure 3. Resilience by excitability heterogeneity** The distance *D* taken from Eq. (22) reflects the degree of the system's resilience: the smaller *D*, the closer are the original and perturbed state and the more resilient is the system towards perturbations. We observe that strong heterogeneity renders the system more resilient to perturbations. For each value of $\sigma_h^2$, 400 values of *D* are plotted to correspond with the same number of simulation trials.</figcaption></figure>

the spectral radius and further limits the number of fixed points. Taken together, these findings suggest that heterogeneity quenches the dynamical repertoire of non-linear networks as their complexity increases to maintain their resilience, suggesting that neural diversity plays a role in the persistence of brain function.

## IV. DISCUSSION

Excitability heterogeneity affects stability of random networks and hence can induce trivialization and enhance resilience. This is shown theoretically in Eq. (5) in gradient systems and in the corresponding application shown in Fig. 3. In non-gradient systems, Equation (11) confirms that an exponential increase in fixed points occurs whenever the May-Wigner instability threshold is met i.e. $\Gamma \ge 1$$^{23-25,28}$. An important corollary of this equation is that the number of equilibria scales with the spectral radius, whose magnitude has been shown to increase with a variety of topological traits, such as increased dimensionality$^{9,18,19}$, strength of connections$^{10-12,14,18,21}$, motifs richness$^{10-12,14,18,21}$, connectance$^{8,9,14,18,21}$ and degree heterogeneity$^{22}$. Interestingly, our results as well as others$^{65}$ indicate that increased intrinsic diversity in cellular excitability represents an exception, contracting the spectral radius and hence promoting stability through a reduction of the number of equilibria. This suggests that variability amongst cells of a network - as opposed to heterogeneity in their connectivity - promotes resilience, by reducing the probability of multistability and bifurcations.

It is interesting to conjecture on the potential role of trivialization in neural circuits' dynamics and function. Intuitively, increasing the number of control parameters by enabling cells with different excitability profiles should enrich a system's phase space – however, our results suggest

the opposite. Numerous experimental$^{2,44-47,49-52}$ and theoretical studies$^{2,53-59}$ have examined the role of cell-to-cell heterogeneity in brain dynamics and function. For instance, several neural population models considered heterogeneous single cell properties$^{90-92}$. Collectively, these indicate that intrinsic diversity may both enhance resilience$^{47,48,52,63,64}$ and optimize coding performance$^{54,56,58,59,62}$. Our complementary results show that variations in cellular excitability represents a stabilizing mechanism, quenching the dynamical repertoire of balanced neural networks. This surprising feature arises notably through the effective linearization of input-output response functions$^{2,23,93}$, promoting graded response to perturbations and subsequent increasingly linear relaxation towards a limited number of stable states.

Multiple mechanisms have been proposed to regularize the dynamics of complex networks beyond heterogeneity. Noise and disorder - ubiquitous in neural systems$^{94}$ - represent important drivers of regularization, where random-like perturbations may paradoxically yield highly predictable dynamics over a wide range of physical$^{67-69}$ and biological$^{70,71}$ systems. This form of complexity collapse$^{95}$ has important implications for neural circuit function. These suggest a trade-off between topological complexity – which has been extensively studied – and intrinsic diversity, balancing coding efficiency and resilience. As such, trivialization resulting from either noise, disorder or heterogeneity, should be considered as a control mechanism, tuning stability towards an optimal dynamical state$^{96-100}$. For instance, destabilizing changes in connectivity during development and learning may be compensated for through adaptive increase in the level of fluctuations, biological diversity or both, to optimize neural circuit performance.

Intrinsic plasticity$^{33-37}$ refers to a non-synaptic form of plasticity through which individual cells' biophysical characteristics and/or excitability change in experience and stimuli-dependent ways, notably through adjustments in ionic conductances$^{47,52,63,101}$. While speculative, our results nonetheless suggest that intrinsic plasticity could play an important role in stabilizing neural dynamics and function, compensating for the notoriously destabilizing effect of synaptic plasticity$^{33,87-89}$.

## DATA AVAILABILITY STATEMENT

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## ACKNOWLEDGMENTS

We thank the National Sciences and Engineering Research Council of Canada (NSERC Grants RGPIN-2017-
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Trivialization may stabilize neural dynamics through development and learning. **A** To model circuit modifications occurring during development and learning, a network of *N* neurons was examined as its connection probability *ρ* increased linearly from 0 to 0.05. **B** In the absence of excitability heterogeneity ($\sigma_h^2 = 0$), the network gradually converges towards instability. As the connection probability increases passed the critical point (vertical dashed line), the network dynamics become unstable. **C** This does not occur if heterogeneity is present in the network ($\sigma_h^2 = 0.001$) and the dynamics remain stable throughout. **D** Spectral disks for the homogeneous ($\sigma_h^2 = 0$; black circle) and heterogeneous excitability ($\sigma_h^2 = 0.001$; blue circle) calculated at intervals of 500 a.u. through Eq. (19). The spectral disks' radius remains below the May-Wigner transition threshold ($\Gamma = 1$; red circle) only in the heterogeneous case. **E** Expected number of equilibria ($\langle N \rangle$) for the homogeneous ($\sigma_h^2 = 0$; black line) and heterogeneous ($\sigma_h^2 = 0.001$; blue line) network as the connection probability increases. In the presence of heterogeneity, the number of equilibria remains small through trivialization, and increases exponentially otherwise. **F** Spectral radius as a function of the variance of the synaptic weight matrix *W*, reflecting destabilizing influence of increase in connectivity and the compensatory effect of heterogeneity (blue lines). One notices that the May-Wigner transition threshold ($\Gamma = 1$; horizontal grey line) is delayed. **G** Associated reduction in the number of equilibria ($\langle N \rangle$) through trivialization at increasing degrees of excitability heterogeneity $\sigma_h^2$ (blue lines). Aside from the connection probability $\rho$ which is increasing, all other parameters are identical to those used in Fig. 2.</figcaption></figure>
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Trivialization may stabilize neural dynamics through development and learning. **A** To model circuit modifications occurring during development and learning, a network of *N* neurons was examined as its connection probability *ρ* increased linearly from 0 to 0.05. **B** In the absence of excitability heterogeneity ($\sigma_h^2 = 0$), the network gradually converges towards instability. As the connection probability increases passed the critical point (vertical dashed line), the network dynamics become unstable. **C** This does not occur if heterogeneity is present in the network ($\sigma_h^2 = 0.001$) and the dynamics remain stable throughout. **D** Spectral disks for the homogeneous ($\sigma_h^2 = 0$; black circle) and heterogeneous excitability ($\sigma_h^2 = 0.001$; blue circle) calculated at intervals of 500 a.u. through Eq. (19). The spectral disks' radius remains below the May-Wigner transition threshold ($\Gamma = 1$; red circle) only in the heterogeneous case. **E** Expected number of equilibria ($\langle N \rangle$) for the homogeneous ($\sigma_h^2 = 0$; black line) and heterogeneous ($\sigma_h^2 = 0.001$; blue line) network as the connection probability increases. In the presence of heterogeneity, the number of equilibria remains small through trivialization, and increases exponentially otherwise. **F** Spectral radius as a function of the variance of the synaptic weight matrix *W*, reflecting destabilizing influence of increase in connectivity and the compensatory effect of heterogeneity (blue lines). One notices that the May-Wigner transition threshold ($\Gamma = 1$; horizontal grey line) is delayed. **G** Associated reduction in the number of equilibria ($\langle N \rangle$) through trivialization at increasing degrees of excitability heterogeneity $\sigma_h^2$ (blue lines). Aside from the connection probability $\rho$ which is increasing, all other parameters are identical to those used in Fig. 2.</figcaption></figure>
INRIA (Action Exploratoire *A/D Drugs* to AH) for support of this research.

## REFERENCES

¹H. Kitano, “Biological robustness,” Nature Reviews Genetics **5**, 826–837 (2004).

²S. Rich, H. M. Chamej, H. Lefebvre, and T. A. Valiante, “Loss of neuronal heterogeneity in epileptogenic human tissue impairs network resilience to sudden changes in synchrony,” Cell Reports **39**, 110863 (2022).

³M. Scheffer, “Critical transitions in nature and society,” in *Critical Transitions in Nature and Society* (Princeton University Press, 2020).

⁴N. Taleb, “Antifragile: Things that gain from disorder,” in *Antifragile: Things That Gain From Disorder* (Random House, 2012).

⁵M. Scheffer, S. R. Carpenter, T. M. Lenton, J. Bascompte, W. Brock, V. Dakos, J. Van de Koppel, I. A. Van de Leemput, S. A. Levin, E. H. Van Nes, et al., “Anticipating critical transitions,” Science **338**, 344–348 (2012).

⁶D. Manheim, “The fragile world hypothesis: Complexity, fragility, and systemic existential risk,” Futures **122**, 102570 (2020).

⁷C. S. Holling, “Resilience and stability of ecological systems,” *Annual review of ecology and systematics* **4**, 1–23 (1973).

⁸R. M. May, “Stability in multispecies community models,” *Mathematical Biosciences* **12**, 59–79 (1971).

⁹R. M. May, “Will a large complex system be stable?” Nature **238**, 413–414 (1972).

¹⁰R. M. May, “Qualitative stability in model ecosystems,” Ecology **54**, 638–641 (1973).

¹¹S. A. Levin, *Some mathematical questions in biology* (American Mathematical Soc., 1974).

¹²R. Levins, “Some demographic and genetic consequences of environmental heterogeneity for biological control,” American Entomologist **15**, 237–240 (1969).

¹³J. Gao, B. Barzel, and A. Barabási, “Universal resilience patterns in complex networks,” Nature **530**, 307–312 (2016).

¹⁴S. Allesina and S. Tang, “Stability criteria for complex ecosystems,” Nature **483**, 205–208 (2012).

¹⁵S. Tang, S. Pawar, and S. Allesina, “Correlation between interaction strengths drives stability in large ecological networks,” Ecology letters **17**, 1094–1100 (2014).

¹⁶L. Fraccascia, I. Giannoccaro, and V. Albino, “Resilience of complex systems: State of the art and directions for future research,” Complexity **2018** (2018).

¹⁷T. Elmqvist, C. Folke, M. Nyström, G. Peterson, J. Bengtsson, B. Walker, and J. Norberg, “Response diversity, ecosystem change, and resilience,” Frontiers in Ecology and the Environment **1**, 488–494 (2003).

¹⁸K. S. McCann, “The diversity–stability debate,” Nature **405**, 228–233 (2000).

¹⁹L. Stone, “The feasibility and stability of large complex biological networks: a random matrix approach,” Scientific Reports **8**, 1–12 (2018).

²⁰J. Ipsen and A. Peterson, “Consequences of dale’s law on the stability-complexity relationship of random neural networks,” Phys. Rev. E **101**, 052412 (2020).

²¹T. Gross, L. Rudolf, S. A. Levin, and U. Dieckmann, “Generalized models reveal stabilizing factors in food webs,” Science **325**, 747–750 (2009).

²²G. Yan, N. Martinez, and Y.-Y. Liu, “Degree heterogeneity and stability of ecological networks,” J. R. Soc. Interface **14**, 20170189 (2017).

²³Y. Fyodorov and B. Khoruzhenko, “Non-linear analogue of the maywigner instability transition,” PNAS **113**, 6827–6832 (2016).

²⁴Y. Fyodorov, “Topology trivialization transition in random non-gradient autonomous odes on a sphere,” Journal of Statistical Mechanics: Theory and Experiment **12**, 124003 (2016).

²⁵Y. Fyodorov and P. Le Doussal, “Journal of statistical physics volume 154, pages466–490 (2014) topology trivialization and large deviations for the minimum in the simplest random optimization,” Journal of Statistical Physics **154**, 466–490 (2014).

²⁶A. R. Ives and S. R. Carpenter, “Stability and diversity of ecosystems,” science **317**, 58–62 (2007).

²⁷P. Landi, H. O. Minoarivelo, Å. Brännström, C. Hui, and U. Dieckmann, “Complexity and stability of adaptive ecological networks: a survey of the theory in community ecology,” in *Systems analysis approach for complex global challenges* (Springer, 2018) pp. 209–248.

²⁸G. Wainrib and J. Touboul, “Topological and dynamical complexity if random neural networks,” Physical Review Letters **110**, 118101 (2013).

²⁹B. Lacroix-A-Chez-Toine and Y. Fyodorov, “Counting equilibria in a random non-gradient dynamics with heterogeneous relaxation rates,” Journal of Physics A: Mathematical and Theoretical **55**, 144001 (2021).

³⁰H. Zeng and J. R. Sanes, “Counting equilibria of large complex systems by instability index,” PNAS **118**, e2023719118 (2021).

³¹R. M. May, *Stability and complexity in model ecosystems* (Princeton University Press, 2001).

³²H. Sompolinsky, A. Crisanti, and H.-J. Sommers, “Chaos in random neural networks,” Physical Review Letters **61**, 259 (1988).

³³G. G. Turrigiano and S. B. Nelson, “Homeostatic plasticity in the developing nervous system,” Nature Reviews Neuroscience **5**, 97–107 (2004).

³⁴W. Zhang and D. J. Linden, “The other side of the engram: Experience-driven changes in neuronal intrinsic excitability,” Nature Rev. Neurosci. **4**, 885 (2003).

³⁵H. Beck and Y. Yaari, “Plasticity of intrinsic neuronal properties in CNS disorders,” Nature Reviews Neuroscience **9**, 357–369 (2008).

³⁶S. J. Tripathy, K. Padmanabhan, R. C. Gerkin, and N. N. Urban, “Intermediate intrinsic diversity enhances neural population coding,” Proceedings of the National Academy of Sciences **110**, 8248–8253 (2013).

³⁷K. Padmanabhan and N. N. Urban, “Intrinsic biophysical diversity decorrelates neuronal firing while increasing information content,” Nature Neuroscience **13**, 1276–1282 (2010).

³⁸M. S. Cembrowski and V. Menon, “Continuous variation within cell types of the nervous system,” Trends in Neurosciences **41**, 337–348 (2018).

³⁹M. S. Cembrowski and N. Spruston, “Heterogeneity within classical cell types is the rule: lessons from hippocampal pyramidal neurons,” Nature Reviews Neuroscience **20**, 193–204 (2019).

⁴⁰F. Scala, D. Kobak, M. Bernabucci, Y. Bernaerts, C. R. Cadwell, J. R. Castro, L. Hartmanis, X. Jiang, S. Laturnus, E. Miranda, et al., “Phenotypic variation of transcriptomic cell types in mouse motor cortex,” Nature **598**, 144–150 (2021).

⁴¹H. Zeng and J. R. Sanes, “Neuronal cell-type classification: challenges, opportunities and the path forward,” Nature Reviews Neuroscience **18**, 530–546 (2017).

⁴²H. Zeng, “What is a cell type and how to define it?” Cell **185**, 2739–2755 (2022).

⁴³B. Tasic, Z. Yao, L. T. Graybuck, K. A. Smith, T. N. Nguyen, D. Bertagnolli, J. Goldy, E. Garren, M. N. Economo, S. Viswanathan, et al., “Shared and distinct transcriptomic cell types across neocortical areas,” Nature **563**, 72–78 (2018).

⁴⁴H. Markram, M. Toledo-Rodriguez, Y. Wang, A. Gupta, G. Silberberg, and C. Wu, “Interneurons of the neocortical inhibitory system,” Nature Reviews Neuroscience **5**, 793–807 (2004).

⁴⁵I. Mody and R. A. Pearce, “Diversity of inhibitory neurotransmission through gabaa receptors,” Trends in neurosciences **27**, 569–575 (2004).

⁴⁶I. Soltesz et al., *Diversity in the neuronal machine: order and variability in interneuronal microcircuits* (Oxford University Press, 2006).
INRIA (Action Exploratoire *A/D Drugs* to AH) for support of this research.

## REFERENCES

¹H. Kitano, "Biological robustness," Nature Reviews Genetics **5**, 826-837 (2004).

²S. Rich, H. M. Chamej, H. Lefebvre, and T. A. Valiante, "Loss of neuronal heterogeneity in epileptogenic human tissue impairs network resilience to sudden changes in synchrony," Cell Reports **39**, 110863 (2022).

³M. Scheffer, "Critical transitions in nature and society," in *Critical Transitions in Nature and Society* (Princeton University Press, 2020).

⁴N. Taleb, "Antifragile: Things that gain from disorder," in *Antifragile: Things That Gain From Disorder* (Random House, 2012).

⁵M. Scheffer, S. R. Carpenter, T. M. Lenton, J. Bascompte, W. Brock, V. Dakos, J. Van de Koppel, I. A. Van de Leemput, S. A. Levin, E. H. Van Nes, et al., "Anticipating critical transitions," Science **338**, 344-348 (2012).

⁶D. Manheim, "The fragile world hypothesis: Complexity, fragility, and systemic existential risk," Futures **122**, 102570 (2020).

⁷C. S. Holling, "Resilience and stability of ecological systems," *Annual review of ecology and systematics* **4**, 1-23 (1973).

⁸R. M. May, "Stability in multispecies community models," *Mathematical Biosciences* **12**, 59-79 (1971).

⁹R. M. May, "Will a large complex system be stable?" Nature **238**, 413-414 (1972).

¹⁰R. M. May, "Qualitative stability in model ecosystems," Ecology **54**, 638-641 (1973).

¹¹S. A. Levin, *Some mathematical questions in biology* (American Mathematical Soc., 1974).

¹²R. Levins, "Some demographic and genetic consequences of environmental heterogeneity for biological control," American Entomologist **15**, 237-240 (1969).

¹³J. Gao, B. Barzel, and A. Barabási, "Universal resilience patterns in complex networks," Nature **530**, 307-312 (2016).

¹⁴S. Allesina and S. Tang, "Stability criteria for complex ecosystems," Nature **483**, 205-208 (2012).

¹⁵S. Tang, S. Pawar, and S. Allesina, "Correlation between interaction strengths drives stability in large ecological networks," Ecology letters **17**, 1094-1100 (2014).

¹⁶L. Fraccascia, I. Giannoccaro, and V. Albino, "Resilience of complex systems: State of the art and directions for future research," Complexity **2018** (2018).

¹⁷T. Elmqvist, C. Folke, M. Nyström, G. Peterson, J. Bengtsson, B. Walker, and J. Norberg, "Response diversity, ecosystem change, and resilience," Frontiers in Ecology and the Environment **1**, 488-494 (2003).

¹⁸K. S. McCann, "The diversity-stability debate," Nature **405**, 228-233 (2000).

¹⁹L. Stone, "The feasibility and stability of large complex biological networks: a random matrix approach," Scientific Reports **8**, 1-12 (2018).

²⁰J. Ipsen and A. Peterson, "Consequences of dale's law on the stability-complexity relationship of random neural networks," Phys. Rev. E **101**, 052412 (2020).

²¹T. Gross, L. Rudolf, S. A. Levin, and U. Dieckmann, "Generalized models reveal stabilizing factors in food webs," Science **325**, 747-750 (2009).

²²G. Yan, N. Martinez, and Y.-Y. Liu, "Degree heterogeneity and stability of ecological networks," J. R. Soc. Interface **14**, 20170189 (2017).

²³Y. Fyodorov and B. Khoruzhenko, "Non-linear analogue of the maywigner instability transition," PNAS **113**, 6827-6832 (2016).

²⁴Y. Fyodorov, "Topology trivialization transition in random non-gradient autonomous odes on a sphere," Journal of Statistical Mechanics: Theory and Experiment **12**, 124003 (2016).

²⁵Y. Fyodorov and P. Le Doussal, "Journal of statistical physics volume 154, pages466-490 (2014) topology trivialization and large deviations for the minimum in the simplest random optimization," Journal of Statistical Physics **154**, 466-490 (2014).

²⁶A. R. Ives and S. R. Carpenter, "Stability and diversity of ecosystems," science **317**, 58-62 (2007).

²⁷P. Landi, H. O. Minoarivelo, Å. Brännström, C. Hui, and U. Dieckmann, "Complexity and stability of adaptive ecological networks: a survey of the theory in community ecology," in *Systems analysis approach for complex global challenges* (Springer, 2018) pp. 209-248.

²⁸G. Wainrib and J. Touboul, "Topological and dynamical complexity if random neural networks," Physical Review Letters **110**, 118101 (2013).

²⁹B. Lacroix-A-Chez-Toine and Y. Fyodorov, "Counting equilibria in a random non-gradient dynamics with heterogeneous relaxation rates," Journal of Physics A: Mathematical and Theoretical **55**, 144001 (2021).

³⁰H. Zeng and J. R. Sanes, "Counting equilibria of large complex systems by instability index," PNAS **118**, e2023719118 (2021).

³¹R. M. May, *Stability and complexity in model ecosystems* (Princeton University Press, 2001).

³²H. Sompolinsky, A. Crisanti, and H.-J. Sommers, "Chaos in random neural networks," Physical Review Letters **61**, 259 (1988).

³³G. G. Turrigiano and S. B. Nelson, "Homeostatic plasticity in the developing nervous system," Nature Reviews Neuroscience **5**, 97-107 (2004).

³⁴W. Zhang and D. J. Linden, "The other side of the engram: Experience-driven changes in neuronal intrinsic excitability," Nature Rev. Neurosci. **4**, 885 (2003).

³⁵H. Beck and Y. Yaari, "Plasticity of intrinsic neuronal properties in CNS disorders," Nature Reviews Neuroscience **9**, 357-369 (2008).

³⁶S. J. Tripathy, K. Padmanabhan, R. C. Gerkin, and N. N. Urban, "Intermediate intrinsic diversity enhances neural population coding," Proceedings of the National Academy of Sciences **110**, 8248-8253 (2013).

³⁷K. Padmanabhan and N. N. Urban, "Intrinsic biophysical diversity decorrelates neuronal firing while increasing information content," Nature Neuroscience **13**, 1276-1282 (2010).

³⁸M. S. Cembrowski and V. Menon, "Continuous variation within cell types of the nervous system," Trends in Neurosciences **41**, 337-348 (2018).

³⁹M. S. Cembrowski and N. Spruston, "Heterogeneity within classical cell types is the rule: lessons from hippocampal pyramidal neurons," Nature Reviews Neuroscience **20**, 193-204 (2019).

⁴⁰F. Scala, D. Kobak, M. Bernabucci, Y. Bernaerts, C. R. Cadwell, J. R. Castro, L. Hartmanis, X. Jiang, S. Laturnus, E. Miranda, et al., "Phenotypic variation of transcriptomic cell types in mouse motor cortex," Nature **598**, 144-150 (2021).

⁴¹H. Zeng and J. R. Sanes, "Neuronal cell-type classification: challenges, opportunities and the path forward," Nature Reviews Neuroscience **18**, 530-546 (2017).

⁴²H. Zeng, "What is a cell type and how to define it?" Cell **185**, 2739-2755 (2022).

⁴³B. Tasic, Z. Yao, L. T. Graybuck, K. A. Smith, T. N. Nguyen, D. Bertagnolli, J. Goldy, E. Garren, M. N. Economo, S. Viswanathan, et al., "Shared and distinct transcriptomic cell types across neocortical areas," Nature **563**, 72-78 (2018).

⁴⁴H. Markram, M. Toledo-Rodriguez, Y. Wang, A. Gupta, G. Silberberg, and C. Wu, "Interneurons of the neocortical inhibitory system," Nature Reviews Neuroscience **5**, 793-807 (2004).

⁴⁵I. Mody and R. A. Pearce, "Diversity of inhibitory neurotransmission through gabaa receptors," Trends in neurosciences **27**, 569-575 (2004).

⁴⁶I. Soltesz et al., *Diversity in the neuronal machine: order and variability in interneuronal microcircuits* (Oxford University Press, 2006).
<footnote n="47">⁴⁷E. Marder and J.-M. Goaillard, "Variability, compensation and homeostasis in neuron and network function," Nature Reviews Neuroscience **7**, 563-574 (2006).</footnote>

<footnote n="48">⁴⁸L. S. Tang, A. L. Taylor, A. Rinberg, and E. Marder, "Robustness of a rhythmic circuit to short-and long-term temperature changes," *Journal of Neuroscience* **32**, 10075-10085 (2012).</footnote>

<footnote n="49">⁴⁹G. S. Tomassy, D. R. Berger, H.-H. Chen, N. Kasthuri, K. J. Hayworth, A. Vercelli, H. S. Seung, J. W. Lichtman, and P. Arlotta, "Distinct profiles of myelin distribution along single axons of pyramidal neurons in the neocortex," *Science* **344**, 319-324 (2014).</footnote>

<footnote n="50">⁵⁰I. Matias, J. Morgado, and F. C. A. Gomes, "Astrocyte heterogeneity: impact to brain aging and disease," *Frontiers in aging neuroscience*, 5 (2019).</footnote>

<footnote n="51">⁵¹H. Moradi Chameh, S. Rich, L. Wang, F.-D. Chen, L. Zhang, P. L. Carlen, S. J. Tripathy, and T. A. Valiante, "Diversity amongst human cortical pyramidal neurons revealed via their sag currents and frequency preferences," *Nature Communications* **12** (2021), 10.1038/s41467-021-22741-9.</footnote>

<footnote n="52">⁵²J.-M. Goaillard and E. Marder, "Ion channel degeneracy, variability, and covariation in neuron and circuit resilience," *Annual review of neuroscience* **44**, 335-357 (2021).</footnote>

<footnote n="53">⁵³C. Börgers and N. Kopell, "Synchronization in networks of excitatory and inhibitory neurons with sparse, random connectivity," *Neural computation* **15**, 509-538 (2003).</footnote>

<footnote n="54">⁵⁴J. Mejias and A. Longtin, "Optimal heterogeneity for coding in spiking neural networks," *Physical Review Letters* **108**, 228102 (2012).</footnote>

<footnote n="55">⁵⁵M. Y. Yim, A. Aertsen, and S. Rotter, "Impact of intrinsic biophysical diversity on the activity of spiking neurons," *Physical Review E* **87**, 032710 (2013).</footnote>

<footnote n="56">⁵⁶J. Mejias and A. Longtin, "Differential effects of excitatory and inhibitory heterogeneity on the gain and asynchronous state of sparse cortical networks," *Front. Comput. Neurosci.* **8**, 107 (2014).</footnote>

<footnote n="57">⁵⁷J. Gjorgjieva, G. Drion, and E. Marder, "Computational implications of biophysical diversity and multiple timescales in neurons and synapses for circuit performance," *Current opinion in neurobiology* **37**, 44-52 (2016).</footnote>

<footnote n="58">⁵⁸R. Gast, S. A. Solla, and A. Kennedy, "Effects of neural heterogeneity on spiking neural network dynamics," *arXiv preprint arXiv:2206.08813* (2022).</footnote>

<footnote n="59">⁵⁹R. Gast, S. A. Solla, and A. Kennedy, "Macroscopic dynamics of neural networks with heterogeneous spiking thresholds," *arXiv preprint arXiv:2209.03501* (2022).</footnote>

<footnote n="60">⁶⁰J. Lengler, F. Jug, and A. Steger, "Reliable neuronal systems: the importance of heterogeneity," *PloS one* **8**, e06094 (2013).</footnote>

<footnote n="61">⁶¹I. Aradi and I. Soltesz, "Modulation of network behaviour by changes in variance in interneuronal properties," *The Journal of physiology* **538**, 227-251 (2002).</footnote>

<footnote n="62">⁶²M. Di Volo and A. Destexhe, "Optimal responsiveness and information flow in networks of heterogeneous neurons," *Sci. Rep.* **11**, 17611 (2021).</footnote>

<footnote n="63">⁶³E. Marder, S. Kedia, and E. Morozova, "New insights from small rhythmic circuits," *Curr Opin Neurobiol.* **76**, 102610 (2022).</footnote>

<footnote n="64">⁶⁴E. Marder, "Neuromodulation of neuronal circuits: back to the future," *Neuron* **76**, 1-11 (2012).</footnote>

<footnote n="65">⁶⁵A. Hutt, S. Rich, T. Valiante, and J. Lefebvre, "Intrinsic neural diversity quenches the dynamic volatility of neural networks," *Proc. Nat. Acad. Sci. USA* **120**, e2218841120 (2023).</footnote>

<footnote n="66">⁶⁶P. S. Skardal, D. Taylor, and J. Sun, "Optimal synchronization of complex networks," *Phys. Rev. Lett.* **113**, 144101 (2014).</footnote>

<footnote n="67">⁶⁷P. S. Skardal and A. Arenas, "Disorder induces explosive synchronization," *Phys. Rev. E* **89**, 062811 (2014).</footnote>

<footnote n="68">⁶⁸G. Hermann and J. Touboul, "Heterogeneous connections induce oscillations in large-scale networks," *Phys. Rev. Lett.* **109**, 018702 (2012).</footnote>

<footnote n="69">⁶⁹F. Scafuti, T. Aoki, and M. di Bernardo, "Heterogeneity induces emergent functional networks for synchronization," *Phys. Rev. E* **91**, 062913 (2015).</footnote>

<footnote n="70">⁷⁰J. D. Touboul, C. Piette, L. Venance, and G. B. Ermentrout, "Noise-induced synchronization and antiresonance in interacting excitable systems: Applications to deep brain stimulation in parkinson's disease," *Phys. Rev. X* **10**, 011073 (2020).</footnote>

<footnote n="71">⁷¹C. Gu, X. Liang, H. Yang, and J. Rohling, "Heterogeneity induces rhythms of weakly coupled circadian neurons," *Sci. Rep.* **6**, 21412 (2016).</footnote>

<footnote n="72">⁷²P. Villegas, J. Ruiz-Franco, J. Hidalgo, and M. A. Munoz, "Intrinsic noise and deviations from criticality in boolean gene-regulatory networks," *Sci. Rep.* **6**, 34743 (2016).</footnote>

<footnote n="73">⁷³K. Rajan and L. Abbott, "Eigenvalue spectra of random matrices for neural networks," *Phys. Rev. Lett.* **97**, 188104 (2006).</footnote>

<footnote n="74">⁷⁴C. van Vreeswijk and H. Sompolinsky, "Chaos in neuronal networks with balanced excitatory and inhibitory activity," *Science* **274**, 1724-1726 (1996).</footnote>

<footnote n="75">⁷⁵V. K. Jirsa and M. Ding, "Will a large complex system with time delays be stable?" *Physical Review Letters* **93**, 070602 (2004).</footnote>

<footnote n="76">⁷⁶V. L. Girko, "Circular law," *Theory of Probability & Its Applications* **29**, 694-706 (1985).</footnote>

<footnote n="77">⁷⁷E. Wigner, "Characteristic vectors bordered matrices with infinite dimensions," *Ann. Math.* **62**, 548-564 (1955).</footnote>

<footnote n="78">⁷⁸J. Ginibre, "Statistical ensembles of complex, quaternion, and real matrices," *Journal of Mathematical Physics* **6**, 440-449 (1965).</footnote>

<footnote n="79">⁷⁹A. Hutt, J. Lefebvre, D. Hight, and H. Kaiser, "Phase coherence induced by additive Gaussian and non-Gaussian noise in excitable networks with application to burst suppression-like brain signals," *Front. Appl. Math. Stat.* **5**, 69 (2020).</footnote>

<footnote n="80">⁸⁰A. Hutt, T. Wahl, N. Voges, J. Hausmann, and J. Lefebvre, "Coherence resonance in random Erdos-Renyi neural networks: mean-field theory," *Front. Appl. Math. Stat.* **7**, 697904 (2021).</footnote>

<footnote n="81">⁸¹J. Lefebvre, A. Hutt, J. Knebel, K. Whittingstall, and M. Murray, "Stimulus statistics shape oscillations in nonlinear recurrent neural networks," *J. Neurosci.* **35**, 2895-2903 (2015).</footnote>

<footnote n="82">⁸²I. Gradshteyn and I. Ryzhik, *Table of Integrals, Series, and Products*, 7th ed. (Academic Press, San Diego, 2007).</footnote>

<footnote n="83">⁸³J. Lefebvre and A. Hutt, "Induced synchronisation by endogenous noise modulation in finite-size random neural networks: a stochastic mean-field study," *HAL hal-04162581 submitted* (2023).</footnote>

<footnote n="84">⁸⁴M. Kac, "On the average number of real roots of a random algebraic equation," *Bulletin of the American Mathematical Society* **49**, 314-320 (1943).</footnote>

<footnote n="85">⁸⁵M. A. Lynch, "Long-term potentiation and memory," *Physiological reviews* **84**, 87-136 (2004).</footnote>

<footnote n="86">⁸⁶J. Aljadeff, M. Stern, and T. Sharpee, "Transition to chaos in random networks with cell-type-specific connectivity," *Phys. Rev. Lett.* **114**, 088101 (2015).</footnote>

<footnote n="87">⁸⁷F. Zenke, G. Hennequin, and W. Gerstner, "Synaptic plasticity in neural networks needs homeostasis with a fast rate detector," *PLoS computational biology* **9**, e1003330 (2013).</footnote>

<footnote n="88">⁸⁸F. Zenke and W. Gerstner, "Hebbian plasticity requires compensatory processes on multiple timescales," *Philosophical Transactions of the Royal Society B: Biological Sciences* **372**, 20160259 (2017).</footnote>

<footnote n="89">⁸⁹W. C. Abraham, "Metaplasticity: tuning synapses and networks for plasticity," *Nature Reviews Neuroscience* **9**, 387-387 (2008).</footnote>

<footnote n="90">⁹⁰D. Goldobin, M. di Volo, and A. Torcini, "Reduction methodology for fluctuation driven population dynamics," *Phys. Rev. Lett.* **127**, 038301 (2021).</footnote>

<footnote n="91">⁹¹R. Stefanescu and V. Jirsa, "A low dimensional description of globally coupled heterogeneous neural networks of excitatory and inhibitory neurons," *PLoS Comput. Biol.* **4**, e1000219 (2008).</footnote>

<footnote n="92">⁹²Y. Zerlaut, B. Telenczuk, C. Deleuze, T. Bal, G. Ouannounou, and A. Destexhe, "Heterogeneous firing rate response of mouse layer v pyramidal neurons in the fluctuation-driven regime," *J. Physiol.* **594**, 3791-808 (2016).</footnote>

<footnote n="93">⁹³S. Rich, A. Hutt, F. K. Skinner, T. A. Valiante, and J. Lefebvre, "Neurostimulation stabilizes spiking neural networks by disrupting seizure-like oscillatory transitions," *Scientific Reports* **10**, 1-17 (2020).</footnote>
<footnote n="47">⁴⁷E. Marder and J.-M. Goaillard, "Variability, compensation and homeostasis in neuron and network function," Nature Reviews Neuroscience **7**, 563-574 (2006).</footnote>

<footnote n="48">⁴⁸L. S. Tang, A. L. Taylor, A. Rinberg, and E. Marder, "Robustness of a rhythmic circuit to short-and long-term temperature changes," Journal of Neuroscience **32**, 10075-10085 (2012).</footnote>

<footnote n="49">⁴⁹G. S. Tomassy, D. R. Berger, H.-H. Chen, N. Kasthuri, K. J. Hayworth, A. Vercelli, H. S. Seung, J. W. Lichtman, and P. Arlotta, "Distinct profiles of myelin distribution along single axons of pyramidal neurons in the neocortex," Science **344**, 319-324 (2014).</footnote>

<footnote n="50">⁵⁰I. Matias, J. Morgado, and F. C. A. Gomes, "Astrocyte heterogeneity: impact to brain aging and disease," Frontiers in aging neuroscience, 59 (2019).</footnote>

<footnote n="51">⁵¹H. Moradi Chameh, S. Rich, L. Wang, F.-D. Chen, L. Zhang, P. L. Carlen, S. J. Tripathy, and T. A. Valiante, "Diversity amongst human cortical pyramidal neurons revealed via their sag currents and frequency preferences," Nature Communications **12** (2021), 10.1038/s41467-021-22741-9.</footnote>

<footnote n="52">⁵²J.-M. Goaillard and E. Marder, "Ion channel degeneracy, variability, and covariation in neuron and circuit resilience," Annual review of neuroscience **44**, 335-357 (2021).</footnote>

<footnote n="53">⁵³C. Börgers and N. Kopell, "Synchronization in networks of excitatory and inhibitory neurons with sparse, random connectivity," Neural computation **15**, 509-538 (2003).</footnote>

<footnote n="54">⁵⁴J. Mejias and A. Longtin, "Optimal heterogeneity for coding in spiking neural networks," Physical Review Letters **108**, 228102 (2012).</footnote>

<footnote n="55">⁵⁵M. Y. Yim, A. Aertsen, and S. Rotter, "Impact of intrinsic biophysical diversity on the activity of spiking neurons," Physical Review E **87**, 032710 (2013).</footnote>

<footnote n="56">⁵⁶J. Mejias and A. Longtin, "Differential effects of excitatory and inhibitory heterogeneity on the gain and asynchronous state of sparse cortical networks," Front. Comput. Neurosci. **8**, 107 (2014).</footnote>

<footnote n="57">⁵⁷J. Gjorgjieva, G. Drion, and E. Marder, "Computational implications of biophysical diversity and multiple timescales in neurons and synapses for circuit performance," Current opinion in neurobiology **37**, 44-52 (2016).</footnote>

<footnote n="58">⁵⁸R. Gast, S. A. Solla, and A. Kennedy, "Effects of neural heterogeneity on spiking neural network dynamics," arXiv preprint arXiv:2206.08813 (2022).</footnote>

<footnote n="59">⁵⁹R. Gast, S. A. Solla, and A. Kennedy, "Macroscopic dynamics of neural networks with heterogeneous spiking thresholds," arXiv preprint arXiv:2209.03501 (2022).</footnote>

<footnote n="60">⁶⁰J. Lengler, F. Jug, and A. Steger, "Reliable neuronal systems: the importance of heterogeneity," PloS one **8**, e06094 (2013).</footnote>

<footnote n="61">⁶¹I. Aradi and I. Soltesz, "Modulation of network behaviour by changes in variance in interneuronal properties," The Journal of physiology **538**, 227-251 (2002).</footnote>

<footnote n="62">⁶²M. Di Volo and A. Destexhe, "Optimal responsiveness and information flow in networks of heterogeneous neurons," Sci. Rep. **11**, 17611 (2021).</footnote>

<footnote n="63">⁶³E. Marder, S. Kedia, and E. Morozova, "New insights from small rhythmic circuits," Curr Opin Neurobiol. **76**, 102610 (2022).</footnote>

<footnote n="64">⁶⁴E. Marder, "Neuromodulation of neuronal circuits: back to the future," Neuron **76**, 1-11 (2012).</footnote>

<footnote n="65">⁶⁵A. Hutt, S. Rich, T. Valiante, and J. Lefebvre, "Intrinsic neural diversity quenches the dynamic volatility of neural networks," Proc. Nat. Acad. Sci. USA **120**, e2218841120 (2023).</footnote>

<footnote n="66">⁶⁶P. S. Skardal, D. Taylor, and J. Sun, "Optimal synchronization of complex networks," Phys. Rev. Lett. **113**, 144101 (2014).</footnote>

<footnote n="67">⁶⁷P. S. Skardal and A. Arenas, "Disorder induces explosive synchronization," Phys. Rev. E **89**, 062811 (2014).</footnote>

<footnote n="68">⁶⁸G. Hermann and J. Touboul, "Heterogeneous connections induce oscillations in large-scale networks," Phys. Rev. Lett. **109**, 018702 (2012).</footnote>

<footnote n="69">⁶⁹F. Scafuti, T. Aoki, and M. di Bernardo, "Heterogeneity induces emergent functional networks for synchronization," Phys. Rev. E **91**, 062913 (2015).</footnote>

<footnote n="70">⁷⁰J. D. Touboul, C. Piette, L. Venance, and G. B. Ermentrout, "Noise-induced synchronization and antiresonance in interacting excitable systems: Applications to deep brain stimulation in parkinson's disease," Phys. Rev. X **10**, 011073 (2020).</footnote>

<footnote n="71">⁷¹C. Gu, X. Liang, H. Yang, and J. Rohling, "Heterogeneity induces rhythms of weakly coupled circadian neurons," Sci. Rep. **6**, 21412 (2016).</footnote>

<footnote n="72">⁷²P. Villegas, J. Ruiz-Franco, J. Hidalgo, and M. A. Munoz, "Intrinsic noise and deviations from criticality in boolean gene-regulatory networks," Sci. Rep. **6**, 34743 (2016).</footnote>

<footnote n="73">⁷³K. Rajan and L. Abbott, "Eigenvalue spectra of random matrices for neural networks," Phys. Rev. Lett. **97**, 188104 (2006).</footnote>

<footnote n="74">⁷⁴C. van Vreeswijk and H. Sompolinsky, "Chaos in neuronal networks with balanced excitatory and inhibitory activity," Science **274**, 1724-1726 (1996).</footnote>

<footnote n="75">⁷⁵V. K. Jirsa and M. Ding, "Will a large complex system with time delays be stable?" Physical Review Letters **93**, 070602 (2004).</footnote>

<footnote n="76">⁷⁶V. L. Girko, "Circular law," Theory of Probability & Its Applications **29**, 694-706 (1985).</footnote>

<footnote n="77">⁷⁷E. Wigner, "Characteristic vectors bordered matrices with infinite dimensions," Ann. Math. **62**, 548-564 (1955).</footnote>

<footnote n="78">⁷⁸J. Ginibre, "Statistical ensembles of complex, quaternion, and real matrices," Journal of Mathematical Physics **6**, 440-449 (1965).</footnote>

<footnote n="79">⁷⁹A. Hutt, J. Lefebvre, D. Hight, and H. Kaiser, "Phase coherence induced by additive Gaussian and non-Gaussian noise in excitable networks with application to burst suppression-like brain signals," Front. Appl. Math. Stat. **5**, 69 (2020).</footnote>

<footnote n="80">⁸⁰A. Hutt, T. Wahl, N. Voges, J. Hausmann, and J. Lefebvre, "Coherence resonance in random Erdos-Renyi neural networks: mean-field theory," Front. Appl. Math. Stat. **7**, 697904 (2021).</footnote>

<footnote n="81">⁸¹J. Lefebvre, A. Hutt, J. Knebel, K. Whittingstall, and M. Murray, "Stimulus statistics shape oscillations in nonlinear recurrent neural networks," J. Neurosci. **35**, 2895-2903 (2015).</footnote>

<footnote n="82">⁸²I. Gradshteyn and I. Ryzhik, Table of Integrals, Series, and Products, 7th ed. (Academic Press, San Diego, 2007).</footnote>

<footnote n="83">⁸³J. Lefebvre and A. Hutt, "Induced synchronisation by endogenous noise modulation in finite-size random neural networks: a stochastic mean-field study," HAL hal-04162581 submitted (2023).</footnote>

<footnote n="84">⁸⁴M. Kac, "On the average number of real roots of a random algebraic equation," Bulletin of the American Mathematical Society **49**, 314-320 (1943).</footnote>

<footnote n="85">⁸⁵M. A. Lynch, "Long-term potentiation and memory," Physiological reviews **84**, 87-136 (2004).</footnote>

<footnote n="86">⁸⁶J. Aljadeff, M. Stern, and T. Sharpee, "Transition to chaos in random networks with cell-type-specific connectivity," Phys. Rev. Lett. **114**, 088101 (2015).</footnote>

<footnote n="87">⁸⁷F. Zenke, G. Hennequin, and W. Gerstner, "Synaptic plasticity in neural networks needs homeostasis with a fast rate detector," PLoS computational biology **9**, e1003330 (2013).</footnote>

<footnote n="88">⁸⁸F. Zenke and W. Gerstner, "Hebbian plasticity requires compensatory processes on multiple timescales," Philosophical Transactions of the Royal Society B: Biological Sciences **372**, 20160259 (2017).</footnote>

<footnote n="89">⁸⁹W. C. Abraham, "Metaplasticity: tuning synapses and networks for plasticity," Nature Reviews Neuroscience **9**, 387-387 (2008).</footnote>

<footnote n="90">⁹⁰D. Goldobin, M. di Volo, and A. Torcini, "Reduction methodology for fluctuation driven population dynamics," Phys. Rev. Lett. **127**, 038301 (2021).</footnote>

<footnote n="91">⁹¹R. Stefanescu and V. Jirsa, "A low dimensional description of globally coupled heterogeneous neural networks of excitatory and inhibitory neurons," PLoS Comput. Biol. **4**, e1000219 (2008).</footnote>

<footnote n="92">⁹²Y. Zerlaut, B. Telenczuk, C. Deleuze, T. Bal, G. Ouannounou, and A. Destexhe, "Heterogeneous firing rate response of mouse layer v pyramidal neurons in the fluctuation-driven regime," J. Physiol. **594**, 3791-808 (2016).</footnote>

<footnote n="93">⁹³S. Rich, A. Hutt, F. K. Skinner, T. A. Valiante, and J. Lefebvre, "Neurostimulation stabilizes spiking neural networks by disrupting seizure-like oscillatory transitions," Scientific Reports **10**, 1-17 (2020).</footnote>
<footnote n="94">⁹⁴A. A. Faisal, L. P. Selen, and D. M. Wolpert, "Noise in the nervous system," Nature Reviews Neuroscience **9**, 292–303 (2008).</footnote>

<footnote n="95">⁹⁵S. Tavakoli and A. Longtin, "Multi-delay complexity collapse," Phys. Rev. Res. **2**, 033485 (2020).</footnote>

<footnote n="96">⁹⁶L. Cocchi, L. L. Gollo, A. Zalesky, and M. Breakspear, "Criticality in the brain: A synthesis of neurobiology, models and cognition," *Progress in neurobiology* **158**, 132–152 (2017).</footnote>

<footnote n="97">⁹⁷Z. Ma, G. G. Turrigiano, R. Wessel, and K. B. Hengen, "Cortical circuit dynamics are homeostatically tuned to criticality in vivo," *Neuron* **104**, 655–664 (2019).</footnote>

<footnote n="98">⁹⁸S.-J. Wang, C. C. Hilgetag, and C. Zhou, "Sustained activity in hierarchical modular neural networks: self-organized criticality and</footnote>

<footnote>oscillations,” Frontiers in computational neuroscience **5**, 30 (2011).</footnote>

<footnote n="99">⁹⁹A. Levina, J. M. Herrmann, and T. Geisel, "Dynamical synapses causing self-organized criticality in neural networks," Nature physics **3**, 857–860 (2007).</footnote>

<footnote n="100">¹⁰⁰S. Ratté, Y. Zhu, K. Y. Lee, and S. A. Prescott, "Criticality and degeneracy in injury-induced changes in primary afferent excitability and the implications for neuropathic pain," *Elife* **3**, e02370 (2014).</footnote>

<footnote n="101">¹⁰¹J. Yang, H. Shakil, S. Ratté, and P. SA, "Minimal requirements for a neuron to coregulate many properties and the implications for ion channel correlations and robustness," *ELife* **1**, e72875 (2022).</footnote>
<footnote n="94">⁹⁴A. A. Faisal, L. P. Selen, and D. M. Wolpert, "Noise in the nervous system," Nature Reviews Neuroscience **9**, 292–303 (2008).</footnote>

<footnote n="95">⁹⁵S. Tavakoli and A. Longtin, "Multi-delay complexity collapse," Phys. Rev. Res. **2**, 033485 (2020).</footnote>

<footnote n="96">⁹⁶L. Cocchi, L. L. Gollo, A. Zalesky, and M. Breakspear, "Criticality in the brain: A synthesis of neurobiology, models and cognition," *Progress in neurobiology* **158**, 132–152 (2017).</footnote>

<footnote n="97">⁹⁷Z. Ma, G. G. Turrigiano, R. Wessel, and K. B. Hengen, "Cortical circuit dynamics are homeostatically tuned to criticality in vivo," *Neuron* **104**, 655–664 (2019).</footnote>

<footnote n="98">⁹⁸S.-J. Wang, C. C. Hilgetag, and C. Zhou, "Sustained activity in hierarchical modular neural networks: self-organized criticality and</footnote>

<footnote>oscillations," *Frontiers in computational neuroscience* **5**, 30 (2011).</footnote>

<footnote n="99">⁹⁹A. Levina, J. M. Herrmann, and T. Geisel, "Dynamical synapses causing self-organized criticality in neural networks," *Nature physics* **3**, 857–860 (2007).</footnote>

<footnote n="100">¹⁰⁰S. Ratté, Y. Zhu, K. Y. Lee, and S. A. Prescott, "Criticality and degeneracy in injury-induced changes in primary afferent excitability and the implications for neuropathic pain," *Elife* **3**, e02370 (2014).</footnote>

<footnote n="101">¹⁰¹J. Yang, H. Shakil, S. Ratté, and P. SA, "Minimal requirements for a neuron to coregulate many properties and the implications for ion channel correlations and robustness," *ELife* **1**, e72875 (2022).</footnote>