Stability and structural constraints of random brain networks with excitatory and inhibitory neural populations

Abstract

The stability of brain networks with randomly connected excitatory and inhibitory neural populations is investigated using a simplified physiological model of brain electrical activity. Neural populations are randomly assigned to be excitatory or inhibitory and the stability of a brain network is determined by the spectrum of the network’s matrix of connection strengths. The probability that a network is stable is determined from its spectral density which is numerically determined and is approximated by a spectral distribution recently derived by Rajan and Abbott. The probability that a brain network is stable is maximum when the total connection strength into a population is approximately zero and is shown to depend on the arrangement of the excitatory and inhibitory connections and the parameters of the network. The maximum excitatory and inhibitory input into a structure allowed by stability occurs when the net input equals zero and, in contrast to networks with randomly distributed excitatory and inhibitory connections, substantially increases as the number of connections increases. Networks with the largest excitatory and inhibitory input allowed by stability have multiple marginally stable modes, are highly responsive and adaptable to external stimuli, have the same total input into each structure with minimal variance in the excitatory and inhibitory connection strengths, and have a wide range of flexible, adaptable, and complex behavior.

Appendices

Appendix A: Brain dynamics model

This appendix summarizes the general brain dynamics model and describes the linearization and transfer to Fourier space that leads to Eq. (1) in Section 2.1. Further details about this model can be found in (Robinson et al. 2002, 2004; Wright et al. 2001).

The cell body potential V _a of neurons in neural population a results when synaptic inputs from afferent neurons in other populations are summed after being filtered and smeared out in time as a result of receptor dynamics and passage through the dendritic tree. It approximately obeys the equation

$$ V_a({\bf r},t) = \sum_b V_{\!ab} ({\bf r},t) $$

(39)

where V _ab is defined by the differential equation

$$ \label{eq:am1} D_{\!ab}V_{\!ab}({\bf r},t) = \nu_{\!ab} \phi_b({\bf r},t-\tau_{\!ab}), $$

(40)

$$ \label{eq:am2} D_{\!ab} = \frac{1}{\alpha_{\!ab} \beta_{\!ab}}\frac{d^2}{dt^2} + \left[\frac{1}{\alpha_{\!ab}} + \frac{1}{\beta_{\!ab}} \right] \frac{d}{dt} + 1, $$

(41)

where 1/β _ab and 1/α _ab are the rise and decay times of the cell-body potential produced by an impulse at a dendritic synapse. The right of Eq. (40) involves contributions from other components ϕ _b , delayed by a time τ _ab owing to anatomical separations. The quantity ν _ab  = N _ab s _ab , N _ab is the mean number of synapses from neurons in component b to neurons in component a and s _ab is the time integrated strength of the response in neurons of type a to a unit signal from neurons of type b.

The mean firing rates Q _a of neurons are nonlinearly related to the mean potentials V _a by Q _a (r,t) = S[V _a (r,t)], where S is a sigmoid function that increases from 0 to Q _max as V _a increases from − ∞ to ∞. We use

$$ \label{eq:am3} S[V_a ({\bf r},t)]=\frac{{\rm Q_{max}}}{1+exp[-{V_a ({\bf r},t)- \theta}/\sigma']}, $$

(42)

where θ is the mean neural firing threshold, relative to resting, and $\sigma' \pi / \sqrt{3}$ is its standard deviation.

The mean firing rate generated by each component gives rise to a neural pulse, whose average over short time scales forms a field ϕ _a (r,t) that propagates at a velocity v _a through axons with a characteristic range r _a . These pulses spread out and dissipate if not regenerated. To a good approximation, this type of propagation obeys a damped-wave equation

$$ \label{eq:am4} \left( \frac{1}{\gamma^{2}_{a}}\frac{\partial^2}{\partial t^2} +\frac{2}{\gamma_{a}} \frac{\partial}{\partial t} +1-r_{a}^{2} \nabla^2 \right) \phi_a ({\bf r},t) = S[V_a ({\bf r},t)], $$

(43)

where γ _a  = v _a / r _a and r _a is the mean range of axons a.

Linearization about an assumed steady state gives

$$ \label{eq:am5} D_a \phi_a ({\bf r},t) = \rho_a V_a({\bf r},t), $$

(44)

$$ \label{eq:am6} D_a = \frac{1}{\gamma_{a}^{2}}\frac{\partial^2}{\partial t^2} +\frac{2}{\gamma_{a}} \frac{\partial}{\partial t}+ 1-r_{a}^{2} \nabla^2, $$

(45)

to first order, where ρ _a  = dQ _a (V _a )/dV _a is the slope of Eq. (42), at the assumed steady state value of V _a . Henceforth, the symbols ϕ _a and V _a denote linear perturbations to the steady state values, since there is no possibility of confusion.

Fourier transforming Eqs. (40)–(45) in time and eliminating V _a yields

$$ \label{eq:am7} D_a(\omega) \phi_a({\bf r},\omega)= \sum_b J_{ab}(\omega) \phi_b({\bf r}, \omega), $$

(46)

where

$$ J_{ab}(\omega)=L_{ab}(\omega)G_{ab} e^{i\omega\tau_{ab}}, $$

(47)

$$ D_a (\omega)=(1-i \omega / \gamma_a )^2-r_a^2 \nabla^2, $$

(48)

$$ L_{ab}(\omega)=\frac{\alpha_{ab}\beta_{ab}}{(\alpha_{ab} -i \omega)(\beta_{ab}-i \omega)}, $$

(49)

$$ G_{ab}=\rho_a \nu_{ab}, $$

(50)

where $[L_{ab}(\omega)]^{-1}$ is the temporal Fourier transform of Eq. (41). Equation (46) describes the linear perturbations of ϕ _a in Fourier space about the assumed steady state. It reduces to Eq. (1) in the text under the assumptions described in Section 2.1.

Appendix B: Spectra and stability of brain networks with randomly distributed inhibitory connections

Previously we studied the spectrum and stability of RCNs (Gray and Robinson 2008b). In this appendix we summarize the results in (Gray and Robinson 2008b) so comparisons can be made with the results in Section 4. RCNs have the same overall gain distribution as RPNs with gain matrix entries distributed with a mean g and variance σ ² given by Eqs. (14) and (15), respectively.

The spectrum of RCNs is accurately described by the results of random matrix theory (Füredi and Komlós 1981; Gray and Robinson 2008a; Sommers et al. 1988) and can be approximated by the superposition of two distributions. The first distribution, called the principal distribution, denoted ρ _rp (x), represents the distribution of an eigenvalue in the spectrum which is real and normally distributed with mean ng and variance σ ²; i.e.,

$$ \label{eq:prineqr} \rho_{rp} (x) = \frac{1}{\sigma\sqrt{2\pi}} {\rm exp}\left[-\left(\frac{x-ng}{2\sigma}\right)\right]. $$

(51)

The second distribution called the bulk distribution, denoted ρ _rb (λ), represents the distribution of the other n − 1 eigenvalues. The bulk distribution approximately equals (with equality as n → ∞)

$$ \label{eq:rsomeq} \rho_{rb}(\lambda) = (\pi n\sigma^2 )^{-1}, $$

(52)

if x ² + y ² ≤ nσ ², and 0 otherwise, where λ = Reλ + iImλ = x + iy. Our work in (Gray and Robinson 2008a) showed that if g ≫ σ ² > 0, the principal and bulk distributions will be clearly distinguishable from each other and the principal distribution will be the distribution of the dominant eigenvalue λ ₁. If g ≈ σ ² then the principal and bulk distributions overlap, while if σ ² ≫ g the bulk distribution will completely overlap the principal distribution and the entire spectral distribution is given by Eq. (52).

Equation (52) shows that at least n − 1 eigenvalues of large random gain matrices will be uniformly distributed within a circle centered on the origin with radius $\sigma\sqrt{n}$. The projection of Eq. (52) onto the real axis is given by

$$ \label{eq:wscl} \rho_{rx} (x)= \int \rho_{rb}(\lambda) dy = \frac{2}{\pi n\sigma^2}\sqrt{n\sigma^2 -x^2}, $$

(53)

if $|x| \leq \sigma\sqrt{n}$ and ρ _rx (0) = 0 otherwise; this is Wigner’s semicircle law (Sommers et al. 1988; Wigner 1967).

The probability P _s that a large random brain network is stable equals the probability that all the eigenvalues in its spectrum satisfy Eq. (9). In (Gray and Robinson 2008a) we showed that the expected probability P _s that a random brain network is stable approximately equals

$$ \begin{array}{rcl} \label{eq:rpos} P_{s}&=& \frac{1}{2}\left[1+{\rm erf}\left(\frac{1-ng}{\sigma\sqrt{2}}\right)\right] \\ \\ &&\times\left[ \frac{1}{2} + \frac{1}{\pi} {\rm sin}^{-1}\left( \frac{1}{\sigma\sqrt{n}}\right) + \frac{\sqrt{n\sigma^2 -1}}{n\pi\sigma^2}\right]^{n-1}. \end{array} $$

(54)

The first term of this product is the probability that principal eigenvalue is less than one while the second term is the probability that all the eigenvalues in the bulk spectrum have real part less than one. Equation (54) shows that P _s →0 as n → ∞ if ng > 1 and nσ ² > 1.

Appendix C: Rajan-Abbott distribution potential

In this appendix we describe the spectral density potential function ψ derived by Rajan and Abbott (Rajan and Abbott 2006) who termed ψ a potential because the spectral density only depends on ψ’s derivatives. To avoid confusion we change the notation used in (Rajan and Abbott 2006) so that it matches what we use here. These results are used in Section 3 to theoretically calculate ρ _b (λ) for random networks with inhibitory populations [Eq. (25)].

The potential derived by Rajan and Abbott is

$$ \label{eq:psi} \psi=\ln \left[ \frac{1+q}{q^{1-p_i}} \right] + \frac{r^2(q\epsilon +1)}{q+1}, $$

(55)

where

$$ \begin{array}{rcl} \label{eq:q} q&=&\frac{(1-\epsilon)r^2 + 2(1-p_i)-1}{2p_i}\\ \\ && + \frac{\sqrt{[(1-\epsilon)r^2-1]^2 + 4(1-p_i)(1-\epsilon)r^2}}{2p_i}. \end{array} $$

(56)

The eigenvalue density inside the circle can be calculated using Eq. (20) via

$$ \label{eq:psi1} \psi'=q'\left[\frac{1+\epsilon r^2}{q+1} - \frac{1-p_i}{q}- \frac{r^2(\epsilon q+1)}{(q+1)^2}\right] +\frac{\epsilon q+1}{q+1}, $$

(57)

and

$$ \begin{array}{rcl} \label{eq:psi2} \psi'' & = & q''\left[\frac{1+\epsilon r^2}{q+1} - \frac{1-p_i}{q} - \frac{r^2(\epsilon q+1)}{(q+1)^2}\right] \\\\& & +2q'\left[\frac{\epsilon}{q+1}-\frac{\epsilon q+1}{(q+1)^2}\right] \\\\ & &+ \left(q'\right)^2\left[\frac{-(1+2\epsilon r^2)}{(q+1)^2}+\frac{1-p_i}{q^2} +\frac{2r^2 (\epsilon q+1)}{(q+1)^3}\right], \end{array} $$

(58)

where, from Eq. (56),

$$ \label{eq:q1} q'=\frac{1\!-\!\epsilon}{2p_i}\left[1\!+\!\frac{(1\!-\!\epsilon)r^2 \!-\!1\!+\! 2(1\!-\!p_i)}{\sqrt{[(1\!-\!\epsilon)r^2\!-\!1]^!+\!4(1\!-\!p_i)(1\!-\!\epsilon)r^2}}\right], $$

(59)

and

$$ \begin{array}{rcl} \label{eq:q2} q'' & = &\frac{(1-\epsilon)^2}{2p_i\sqrt{\{(1-\epsilon)r^2-1\}^2+4f(1-\epsilon)r^2}} \\ \\ && \times \left[1-\frac{\{(1-\epsilon)r^2 -1+2(1-p_i)\}^2}{[(1-\epsilon)r^2-1]^2+4(1-p_i) (1-\epsilon)r^2}\right]. \end{array} $$

(60)