Stability constraints on large-scale structural brain networks

Introduction

The brain is possibly the most complicated example of a system of interacting dynamical units whose activity self-organizes to produce complex global behavior. The human brain performs cognitive functions through the transmission of action potentials within a vast structurally dynamic network consisting of approximately 10¹¹ neurons and up to 10¹⁵ synaptic interconnections (Kandel et al., 2000; Koch, 2004; Sporns et al., 2005). The aggregate of all neural firings within this network results in large-scale coherent electrical activity and the performance of high-level cognitive functions. Understanding the structure and physiology of the brain thus gives insight into its overall behavior.

At large scales the excitatory and inhibitory neurons in the brain are organized into a complex large-scale network of distinct anatomical and functional structures (Sporns et al., 2004, 2005; Bullmore and Sporns, 2009; van den Heuvel and Sporns, 2011) We can represent this structure as a complex network—the structure of which has been studied extensively in recent times with a number of experimental cortical connection networks determined for the cat and the macaque monkey (Felleman and van Essen, 1991; Scannell et al., 1995; Jouve et al., 1998; Hilgetag et al., 2000a,b; Sporns, 2011). These networks have a modular hierarchical structure with the small-world properties of high local clustering and short path length between structures (Hilgetag et al., 2000a,b; Sporns et al., 2000, 2004, 2005; Young, 2000; Sporns and Zwi, 2004; Bassett and Bullmore, 2006; Bullmore and Sporns, 2009, 2012).

Reasons for why the brain has evolved this particular large-scale structure are currently unknown. A number of investigations have concentrated on the effect of physical constraints on brain structure. Such constraints include brain volume, wiring length, and energy consumption or metabolic demands (Laughlin et al., 1998; Attwell and Laughlin, 2001; Lennie, 2003). Other studies have looked at functional constraints such as minimizing the conduction delay or processing steps for a signal to travel from one neuron to another (Wen and Chklovskii, 2005). Alternatively, the dynamics of the brain's electrical activity may constrain the brain's structure. If the physiology and structural characteristics of the brain produce adverse electrical activity resulting in seizures, tremors, or other neurological disorders then it is likely the structural characteristics of the brain will be constrained to limit these disorders.

One of the most important dynamical properties of complex systems such as the brain is stability. It has been associated with pattern formation (Turing, 1952; Murray, 2002), synchronized activity (Kuramoto, 1984; Pecora and Carroll, 1998; Jirsa and Ding, 2004; Acebrón et al., 2005; Feng et al., 2006), the complexity and diversity of ecosystems (May, 1972, 1974; Hogg et al., 1989; McCann, 2000; Allesina and Tang, 2012), the functioning of biological systems (Murray, 2002; Taverna and Goldstein, 2002; Steuer, 2007), and the generation of coherent self-organized behavior. Stability is also an important aspect of the design and control of advanced technological systems (Bechhoefer, 2005).

A common approach to studying the large-scale dynamics of the brain's electrical activity is to use a continuum mean-field approximation for neural activity. This approach has been extensively studied over the past 30 years producing numerous models for the electrical activity within the brain (Wilson and Cowan, 1973; da Silva et al., 1974; Nunez, 1974, 1995; Freeman, 1975; Steriade et al., 1990; Jirsa and Haken, 1996; Wright and Liley, 1996; Robinson et al., 1997, 2003a; Wright et al., 2001; Robinson, 2005). This work has been reviewed recently (Deco et al., 2008; Bressloff, 2011). These models have been used extensively to perform stability analysis and understand pattern formation, oscillations and waves in the brain's electrical activity (Deco et al., 2008; Bressloff, 2011).

Measurements of brain activity suggest that the brain operates close to marginal stability, permitting a wide range of flexible, adaptable, and complex behavior (Stam et al., 1999; Robinson et al., 2001b; Breakspear, 2002; Breakspear et al., 2003). Physiological modeling also suggests that linear instabilities in the brain's electrical activity correspond to neurological disorders, such as epilepsy (Robinson et al., 1998, 2002; Breakspear et al., 2006; Kim and Robinson, 2007; Deco et al., 2008). It is therefore possible that stability is a dynamical property that imposes constraints on the brain's physiology and structure.

In previous work we have used a simplified version of the Robinson, Rennie, Wright (RRW) physiologically-based continuum (mean-field) model (Robinson et al., 1997, 1998, 2003a; Wright et al., 2001; Rennie et al., 2002; Robinson, 2003, 2005) to study the dynamics of the electrical activity in large-scale structural brain networks. Based on the hypothesis that stability is a dynamical constraint on the structure and physiology of structural brain networks we have investigated the effect of stability of large-scale structural brain networks (Gray and Robinson, 2006, 2008, 2009b,a; Robinson et al., 2009).

Our previous work ignored the dendritic time constants of neurons and the propagation time delays for signals to travel between neural populations. However, time delays due to axonal propagation affect stability and the possible physiology of neuronal networks (Atay and Hutt, 2004; Jirsa and Ding, 2004; Coombes, 2005; Coombes et al., 2007; Qubbaj and Jirsa, 2007, 2009; Venkov et al., 2007; Jirsa, 2009). This previous work has generally used integro-differential neural field equations with connectivity within a neural mass described by homogeneous or heterogeneous kernels. Our approach here is to focus on the temporal dynamics of the overall electrical activity of arbitrarily connected large-scale structural brain networks, ignoring the spatial spread and propagation of electrical activity within individual neuronal populations.

In this study, we increase the physiological realism of our structural brain network model by allowing propagation time delays and non-zero dendritic time constants. After reviewing the stability of structural brain networks, we aim to investigate how these physiological features might affect the dynamics and stability of networks of neuronal populations where the connection patterns between populations are arbitrary—ignoring the spatial and geometric placement of the populations and simply focusing on which populations are inter-connected.

We also investigate the dispersion frequencies of marginally stable modes of electrical activity using plausible physiological parameters. The incorporation of non-zero dendritic time constants generalizes the work in (Jirsa and Ding, 2004) which, by including time delays, extended May's original analysis on the stability of complex systems (May, 1972, 1974).

Methods

A structural brain network of n neural populations is represented by a directed graph N whose vertices and edges represent specific neural populations and inter-population connections, respectively. Neural populations within a network are collections of neurons with an assumed effective range and of sufficient number for a mean-field approximation to be valid. For example, a neural population can represent all the neurons in a distinct region or nuclei of the brain (e.g., cortical area, thalamus), a particular neuron type (e.g., interneuron, pyramidal cell), or a particular neurotransmitter type (e.g., glutamate, GABA, dopamine). Neurons in one population do not have to be separated geometrically or physically within the brain and can be intermixed with the neurons of another population (e.g., excitatory and inhibitory neurons in the cortex).

The structure of N is represented by a connection matrix C(N) = [C_ab]; where C_ab = 1 if there is a connection from population b to population a, C_ab = 0 otherwise. If C_ab = C_ba for all a and b, the network is symmetric; otherwise it is asymmetric. Self-connections in structural brain networks correspond to non-zero diagonal entries in C(N). The connection matrix simply records whether one neuronal population sends neural signals to another neuronal population. Properties of connections are not included in C(N).

Physiologically-Based Structural Brain Network Dynamics

In this section we outline the physiological model used to describe the dynamics of a brain network. If a neural population contains a sufficient number of neurons a continuum approximation can be used, whereby the properties of population neurons are averaged over. This approximation is valid for length scales greater than a few tenths of a millimeter and is thus suitable for investigating the dynamics of large-scale structural brain networks.

The continuum approximation allows the use of a previously developed model, the RRW model, for the brain's electrical activity (Robinson et al., 1997; Wright et al., 2001; Rennie et al., 2002; Robinson et al., 2003b, 2004; Robinson, 2005). This continuum model incorporates and describes three features of neural dynamics: (1) the synapto-dendritic dynamics resulting in the cell body potential; (2) from the mean cell body potential an average firing rate is determined via a non-linear sigmoid function; and (3) the population firing rate generates a neural pulse forming a field ϕ(t) that propagates along the populations outgoing connections. The field within a population is temporally described using a damped wave equation. Implicitly the neurons in each population are assumed to have an effective range which gives a rate at which spikes reach axonal terminals and cease existance.

This model has been extensively used to model the corticothalamic system with the linear version having been shown to produce excellent agreement with EEG spectra, ERP, and other neurophysical phenomena (Robinson et al., 1997, 2001a,b; Rennie et al., 2002; Robinson, 2003). To apply this continuum model to brain networks we previously used a number of simplifying assumptions. In particular, we assumed that all neural populations have instantaneous dendritic response times and there is no time delay for a signal to be sent from one population to the other (Gray and Robinson, 2006, 2008, 2009b,a; Robinson et al., 2009). For this study we relax some of these assumptions.

Firstly we assume time delays τ for a signal to be sent from one population to another are equal. Secondly, we assume each population in a network has the same dendritic decay rate α and rise rate β. The values of 1/α and 1/β equal the dendritic decay and rise time constants, respectively; instantaneous rise and decay times imply 1/α = 1/β = 0. These assumptions are unrealistic for real structural brain networks but improve our previous analysis and allow us to analytically determine stability. Though for some structural brain networks these assumptions may be good approximations of the networks physiology. Another weakness of these assumptions is self-connections have the same time delay as connections between distinct populations—if τ ≠ 0 then self-connections in a network also involve a delay. Generally, self-connections represent interconnections within a neural population and would be expected to have zero time delay. However, for cortical networks self-connections can be used to represent feedback from underlying structures such as the thalamus. Time delayed self-connections would be appropriate for this type of feedback. We will generalize the assumption of equal time delays in future work.

The neurophysics and neurophysiology incorporated into the general RRW model and the equations for the linear perturbations of the neural field ϕ_a, for each neural population a, are described and derived in detail elsewhere (Robinson et al., 1997; Wright et al., 2001; Rennie et al., 2002; Robinson et al., 2003b, 2004; Robinson, 2005). This study uses the notation and equations derived in (Robinson, 2005). Under the assumptions used here the RRW equations describing linear perturbations of the neural field ϕ_a of population a about the assumed steady state in Fourier space reduce to

$\begin{array}{ll}{\left(1-i\omega /\gamma \right)}^2{\varphi }_a(\omega )={\displaystyle \sum _{b}L}(\omega )G_{ab}{e}^{i\omega \tau }{\varphi }_b(\omega ),\hfill & (1)\hfill \\ =L(\omega )e^{i\omega \tau }{\displaystyle \sum _b{G}_{ab}}{\varphi }_b(\omega ),\hfill & (2)\hfill \end{array}$

where ω is the angular frequency and

$\begin{array}{cc}L(\omega )=\frac{\alpha \beta }{(\alpha -i\omega )(\beta -i\omega )}=\frac{1}{(1-i\omega /\alpha )(1-i\omega /\beta )}.& (3)\end{array}$

The gain G_ab is a dimensionless quantity describing the effect of changes in the firing rate of neurons in population b on the neurons of population a. Physiologically, G_ab is the number of extra action potentials produced in a per extra action potential incident from b. Hence, G_ab is a measure of how sensitive and responsive a is to changes in b's activity. In the general RRW model γ is a damping rate equal to the velocity of the ϕ's propagation within a neural population divided by the characteristic range of the axons that carry it. In the spatially uniform case used here, γ represents a temporal damping rate.

Letting G = [G_ab] be the matrix of gains and setting

$\begin{array}{cc}\begin{array}{c}D(\omega )={\left[L\left(\omega \right)\right]}^{-1}\left(1-i\omega /\gamma \right)e^{-i\omega \tau },\\ =(1-i\omega /\alpha )(1-i\omega /\beta ){\left(1-i\omega /\gamma \right)}^2{e}^{-i\omega \tau },\end{array}& \begin{array}{c}(4)\\ (5)\end{array}\end{array}$

which is a complex analytic function. Equation (1) can be written in matrix form as

$\begin{array}{cc}\Delta (\omega )\Phi (\omega )=G\Phi (\omega ),& (6)\end{array}$

where Φ is a column vector of the ϕ_a and Δ(ω) = D(ω)I, where I is the identity matrix. Setting A = G − Δ, Equation (6) can be simplified to

$\begin{array}{cc}A(\omega )\Phi (\omega )=0.& (7)\end{array}$

The linear stability of a network is then determined by the solutions ω of the dispersion relation,

$\begin{array}{cc}\det [A(\omega )]=0.& (8)\end{array}$

The gain matrix G = [G_ab] encodes all of the information in C, since G_ab ≠ 0 implies C_ab ≠ 0, as well as the strength of connections between populations. No assumptions (such as homogeneity or isotropy) are made for the characteristics of a connection and any attenuation or phase shifting of an incoming signal due to time delays are reflected in the exponential term of Equation (1). However, the model implicitly assumes an effective range for neurons within a population. If G_ab > 0 then the connection is excitatory and if G_ab < 0 the connection is inhibitory. Note that Equation (5) shows that if the values of α and β are exchanged, the brain network has the same dynamics and stability.

Realistic Parameter Values for Large-Scale Structural Brain Networks

Physiologically plausible parameter values for γ, α, β, and τ are shown in Table 1. These values are based on the parameters used in the corticothalamic model (Robinson et al., 2003a, 2004), with the specific values taken from (Robinson et al., 2004). The nominal values in Table 1 are the default model parameters used to illustrate our results under our assumptions α and β are the same value for all populations and τ is the same for all connections. The values for γ are based on the cortical excitatory neurons which form the long range connections within the cortex. Inhibitory inter-neurons in the cortex are short range (Nunez, 1995) and therefore have γ ≈ ∞. Under the assumptions used here all neural populations are given the same γ value.

TABLE 1

Table 1. Physiologically plausible ranges and nominal values of parameters.

In real structural brain networks dendritic time constants and propagation time delays may vary. The spatial distribution and physical separation of structures within the brain will lead to distinct time delays. The values of τ in Table 1 are physiologically plausible values for the time delays based on corticothalamic modeling. The nominal value is a realistic value for the average delay between the large-scale neural populations (or areas) in the cerebral cortex.

Randomly Connected Large-Scale Structural Brain Networks

To illustrate our results we investigate randomly connected structural brain networks where neural populations are connected randomly with probability p. The size n and probability of connection we use is based on experimentally determined cortical connection networks for animals. These have been analyzed with graph-theoretical methods and all of these networks have less than 100 neural populations with a connection density (percentage of existing connections out all possible connections) of 20–40% (Felleman and van Essen, 1991; Scannell et al., 1995; Hilgetag et al., 2000a; Sporns et al., 2000, 2004; Bullmore and Sporns, 2009; Rubinov and Sporns, 2010; Sporns, 2011). We use random networks with n = 50 and p = 0.5 to illustrate our results, allowing comparisons with real cortical networks. These values of p ensure the networks are strongly connected (Bollobás, 1985) and all populations have at least one input and one output with high probability; i.e., there are no sources or sinks of electrical activity.

The specific random networks we investigate are the same random networks we have previously investigated (Gray and Robinson, 2006, 2008, 2009a,b). These networks consist of excitatory and inhibitory connections. The probability that a connection is inhibitory is given by p_i and such a connection has a negative gain. Excitatory gains are given values from a normal distribution with a mean μ_e > 0 and variance σ²_e. Similarly, inhibitory connections have a gain sampled from a normal distribution with μ_i < 0 and variance σ²_i. In terms of the gain matrix G all positive entries are sampled from (μ_e, σ²_e) and all negative entries are sampled from (μ_i, σ²_i).

Based on these parameters we investigated the stability of three types of networks: random networks with fixed excitatory gains (RENs), random connection networks (RCNs) with excitatory and inhibitory connections distributed randomly within the network (Gray and Robinson, 2009b), and random population networks (RPNs) (Gray and Robinson, 2009a). RPNs represent random networks with excitatory and inhibitory populations of neurons, this implies the outgoing connections of a given populations are all excitatory or all inhibitory. The gain matrix of RPNs consists of columns with either all entries ≤ 0 or all entries ≥ 0.

We determine the dispersion solutions for these structural brain networks numerically using a FORTRAN program called CROOT (Botten et al., 1983). This program finds dispersion solutions by implementing a recursive algorithm that employs Cauchy's integral formula (Mitrinović and Keckić, 1984) within a specified annulus or disk.

Results

Our results describe the stability of structural brain networks by determining the criteria for a network to stable—starting from simple excitatory networks and then adding time delays and dendritic time constants. For our network model we show that stability is determined by the eigenvalues of the gain matrix with stability constraining the eigenvalues to a specific zone in the complex plain. The first subsections translate the results from previous work into the current context. In particular, when time delays are included we produce a similar tear-drop shaped stability zone found by others (Jirsa and Ding, 2004; Feng et al., 2006; Qubbaj and Jirsa, 2007; Jirsa, 2009; Qubbaj and Jirsa, 2009). However, we show the addition of dendritic time constants modifies the shape of the stability zone. Finally, we use our results to assess how stability constrains the physiology of randomly connected networks with excitatory and inhibitory connections.

Stability of Structural Brain Networks

The solutions ω of the dispersion relation Equation (8) determine the linear stability of a network. Setting λ = D(ω), the dispersion relation is

$\begin{array}{cc}\det \left(G-\lambda I\right)=0.& (9)\end{array}$

Therefore, network stability is determined by the spectrum of G, which we denote Sp(G). All the dispersion solutions ω of the network can be obtained by solving

$\begin{array}{cc}\lambda -D(\omega )=0& (10)\end{array}$

for each λ in Sp(G). If all the λ in Sp(G) have corresponding ω [given by Equation (10)] with Imω < 0 then the network is stable. However, if there exists one λ which has a corresponding dispersion solution with Imω ≥ 0 then the network is unstable. The set of dispersion solutions of a brain network is termed the dispersion spectrum. Taking the complex conjugate of Equations (4, 10) show that if ω₁ = Reω + iImω = ω_r + iω_i is solution for λ then ω₂ = −ω_r + iω_i is a solution for the complex conjugate $\overline{\lambda }$ of λ. Therefore, since both λ and $\overline{\lambda }$ are in Sp(G), the dispersion spectrum is symmetric about the real axis.

Solving Equation (10) for ω is equivalent to solving λ −(ϖ) = 0 for ϖ where

and ϖ = ω/γ is a dimensionless frequency parameter. From the ϖ solutions the dispersion solutions for the network are ω = γϖ. Since D and are analytic, the dispersion spectrum can be obtained by numerically solving Equation (10) for each λ.

Boundary between unstable and stable states

As the stability of a network is determined by Sp(G) we are interested in the zone of the complex plane where the all the eigenvalues of G must lie for the network to be stable. If a dispersion solution has Imω = 0 (i.e., ω is real and marginally stable) then the λ corresponding to ω lies on the critical boundary between the unstable and stable zones in the complex plane. Therefore, the stability boundary is given by D(ω) for real ω ranging from −∞ to ∞. To describe the stability boundary we use (ϖ) with real ϖ. This function traces out a continuous curve in the complex plane as ϖ ranges from −∞ to ∞. Points on this curve represent λ in Sp(G) with real dispersion solutions ω = γϖ. The stability boundary of a network is a segment of this curve since, in general, (ϖ) can intersect itself to form loops.

Since is a re-parameterization of D using ϖ = ω/γ, a brain network with model parameters γ = γ′, α = α′, β = β′, and τ = τ' has the same stability boundary as a network with γ = 1, α = α′/γ′, β = β′/γ′, and τ = γ′τ. However, the dispersion spectra of these networks will differ. If an initially stable network becomes unstable due to changes in its connection gains, its spectrum initially lies in the stability zone before at least one eigenvalue moves across the stability boundary from the stable to unstable zones. The λ that crosses the boundary is an instability with frequency ω/2π = γϖ/2π, where λ = (ϖ) on the boundary.

We now describe some general properties of the stability boundary for structural brain networks with time delays and non-zero dendritic time constants. In the following sections we investigate particular cases. Firstly, ≠ 0 for all real −∞ < ϖ < ∞, thus λ = 0 does not lie on the stability boundary. If λ = 0 is substituted into Equation (11) then

$\begin{array}{cc}{(1-i\varpi )}^2=0,{(1-i\varpi \gamma /\alpha )}^2=0,\text{or}{(1-i\varpi \gamma /\beta )}^2=0,& (12)\end{array}$

for complex ϖ. The only solutions to Equation (12) are ϖ = −i, ϖ = −iα/γ, or ϖ = −iβ/γ giving the dispersion solution ω = −iγ, ω = −iα, and ω = −iβ. These solutions all lie in the lower half plane for α, β, γ > 0 and hence λ = 0 is a stable eigenvalue. If ω = ϖ = 0 then (ϖ) = 1, since L(0) = 1, and λ = 1 lies on the stability boundary. This implies that if λ = 1 is in Sp(G) then the network has a zero frequency marginally stable dispersion solution. Furthermore, if we consider

$\begin{array}{ll}\left|\lambda \right|=|L(\varpi ){|}^{-1}|e^{-i\varpi \gamma \tau }\left|\right|1-i\varpi {|}^2,\hfill & (13)\hfill \\ =|1-i\varpi \gamma /\alpha ||1-i\varpi \gamma /\beta |\left|e^{-i\varpi \gamma \tau }\right||1-i\varpi {|}^2,\hfill & (14)\hfill \end{array}$

then |λ| < 1 implies at least one of the factors on the right of Equation (14) is less than 1. Since α, β, γ, and τ are all positive, this condition can only be satisfied if Imω < 0. Therefore, the unit disk is always contained in the stability zone and if all the eigenvalues of a gain matrix lie in the unit disk the network is stable independent of α, β, γ, and τ.

Stability of excitatory networks

If the all the connections in a structural brain network are excitatory then G_ab ≥ 0 and G is a non-negative matrix. The Perron–Frobenius theorem (Horn and Johnson, 1985; Cvetković et al., 1995) then implies that G has a real eigenvalue λ_p such that |λ_i| ≤ λ_p for all λ_i in Sp(G). Therefore, an excitatory brain network is stable if and only if λ_p < 1; i.e., all the eigenvalues are in the unit disk. This stability criteria follows from the discussion at the end of the previous section and the fact that if λ_p = 1 then ω = 0 is a solution to Equation (11). Since 0 lies on the stability boundary this implies if λ_p ≥ 1 then Imω ≥ 0. This means the stability of an excitatory brain network is independent of α, β, γ, and τ (as described in the previous section). In general, structural brain networks have inhibitory connections and the Perron–Frobenius theorem does not apply. This implies the presence of inhibitory connections allows the stability zone to extend beyond the unit disk.

Impact of Time Delays and Dendritic Time Constants on Network Stability

Stability of networks with no time delays and instantaneous rise and decay times

For networks with no time delays and instantaneous dendritic rise and decay times (i.e., 1/α = 1/β = τ = 0)

and hence, for each λ in Sp(G) there are two dispersion solutions given by

$\begin{array}{cc}\omega =\gamma \varpi =-\gamma (i\pm i\sqrt{\lambda }).& (16)\end{array}$

Taking the imaginary part of Equation (16) we obtain

$\begin{array}{cc}\text{Im}\omega =-\gamma (1\pm \text{Re}\sqrt{\lambda }).& (17)\end{array}$

The stability condition Imω < 0 implies that −γ(1 ± Re $\sqrt{\lambda }$) or Re$\sqrt{\lambda }$ ≤ 1 since γ > 0. Thus the stability of a network with 1/α = 1/β = τ = 0 is independent of γ and all λ must satisfy Re$\sqrt{\lambda }$ ≤ 1 or alternatively Reλ + |λ| ≤ 2 [since (Re$\sqrt{\lambda }$)² equals (Reλ + |λ|)/2, with equality corresponding to the stability boundary].

If λ_r = Reλ and λ_i = Imλ then the stability zone is a parabolic zone in the complex plane given by

$\begin{array}{cc}{\lambda }_i^2\le 4-4{\lambda }_r& (18)\end{array}$

with equality giving the stability boundary. The axis of the parabolic boundary is along the real axis with a turning point at (λ_r, λ_i) = (1, 0) and imaginary axis intercepts at λ_i = ± 2. This stability zone is the light gray zone in Figure 1 and is the stability region described previously in (Gray and Robinson, 2008, 2009a,b).

FIGURE 1

Figure 1. Stability zone for a brain network with 1/α = 1/β = τ = 0. The gray zone is where all the eigenvalues λ of the brain network must lie for the network to be stable. The dark region within this zone is the unit disk.

Stability of networks with equal time delays

We now consider the effect of time delays on the stability of networks with 1/α = 1/β = 0. In this section we determine how the addition of a time delay to structural brain networks modifies the parabolic stability zone described by Equation (18) and describe the characteristics of the stability boundary. As noted previously, under our assumptions any self-connection in a brain network has the same time delay as connections between populations.

In this case L(ω) = 1, τ ≠ 0, and

Due to the exponential in Equation (19), λ − (ϖ) = 0 has an infinite number of solutions for each eigenvalue of the gain matrix. If λ = Reλ + iImλ = λ_r + iλ_i then λ = (ϖ) (for real ϖ) implies

and

these equations can be combined giving

$\begin{array}{cc}{\lambda }_i=-\tan (\varpi \gamma \tau )[{\lambda }_r+2\varpi \text{cosec}(\varpi \gamma \tau )],& (22)\end{array}$

which gives the stability boundary in the complex plane. For ϖ ≥ 0, Equation (22) describes a spiral curve in the complex plane starting at the point (λ_r, λ_i) = (1, 0), centered on the origin, and spiraling in a clockwise direction as ϖ increases, as seen in Figure 2A. For ϖ < 0, Equation (22) describes a similar counterclockwise spiral curve, corresponding to a reflection of the curve in Figure 2A about the imaginary axis.

FIGURE 2

Figure 2. Stability boundaries and stability zones for networks with 1/α = 1/β = 0. The curve (ϖ) for ϖ ≥ 0 (A) and corresponding stability boundary (B) for networks with 1/α = 1/β = 0 and γτ = 3.0. The shaded zone in (B) is the stability zone for a network and the dark region within this zone is the unit disk. The dot-dashed line is the boundary for 1/α = 1/β = τ = 0. (C) Stability boundaries for brain networks with 1/α = 1/β = 0 and γτ = 0.1, 0.5, 1.0, and 5.0. The values of γτ are written slightly above and to the left of the corresponding boundary. The dark gray region within this zone is the unit disk.

Using Equations (20–22) we now describe the characteristics of the stability boundary and the resulting stability zone. The (ϖ) curve crosses the real axis when λ_i = 0. Substituting λ_i = 0 into Equation (21) and (22) gives

$\begin{array}{cc}({\varpi }^2-1)\text{sin}(\varpi \gamma \tau )=2\varpi \text{cos}(\varpi \gamma \tau ),& (23)\end{array}$

and

$\begin{array}{cc}{\lambda }_r=-2\varpi \text{cosec}(\varpi \gamma \tau ),& (24)\end{array}$

respectively. The values of ϖ for which λ_i = 0 can be obtained by solving Equation (23) numerically. Note that Equation (24) is only valid if ϖγτ ≠ ±mπ/2 for integers m. Due to the periodicity of the sine and cosine functions there is an infinite number of ϖ that satisfy these equations.

Thus the stability boundary consists of two spiral curves produced by bending the arms of the parabola described by Equation (18) inwards. These spiral curves intersect an infinite number of times on the real axis enclosing larger and larger regions of the complex plane as |ϖγτ| increases. The intersection of all these enclosed regions, corresponding to the innermost zone, represents the stability zone for the network. The reason eigenvalues outside this innermost zone are instabilities is because they correspond to eigenvalues outside the stability zone described by Equation (18) when transformed through multiplication by e^iϖγτ (which removes the effect of the time delay).

We define the smallest ω = γϖ > 0 that gives λ_i = 0 the critical ω value. This critical value is denoted by ω_c with the corresponding ϖ and λ_r denoted ϖ_c and λ^c_r, respectively. The stability zone is defined by

for 0 ≤ ϖ ≤ ϖ_c. This zone has a boundary defined by Equation (19) for −ϖ_c ≤ ϖ ≤ ϖ_c. The real axis intercepts of the stability boundary are given by λ_r = 1 and λ_r = λ^c_r = −2ϖ_ccosec(ϖ_cγτ). As seen in Figure 2B the stability zone has a teardrop shape containing the unit disc. Note that Equation (25) also defines the stability zone in Figure 1 with ϖ_c = ∞. As γτ → 0 then λ_r → (1 − ϖ²) in Equation (20), λ_i → −2ϖ in Equation (21), and the stability boundary converges to the parabola λ²_i = 4 − 4λ_r, as expected.

The effect of increasing τ on the stability boundary is shown in Figure 2C. As γτ increases, the stability boundary converges to the unit circle with the stability zone converging to the unit disk, shaded dark gray. The values of ϖ_c in Figure 2C for the four values of γτ shown are 3.0 (γτ = 0.1), 1.92 (γτ = 0.5), 1.31 (γτ = 1.0), and 0.46 (γτ = 5.0), respectively. The corresponding λ^c_r are −9.4 (intersection not seen), −4.6, −2.7, and −1.2. For γ = 100 s⁻¹ and τ = 0.01 s (the nominal values in Table 1) ω_c/2π ≈ 30 Hz.

Overall, the presence of time delays bends the parabola described by Equation (18) inward, forming a teardrop-shaped stability zone containing the unit disk. As τ → ∞ the stability boundary wraps around the unit circle an infinite number of times and the stability zone converges to the unit disk, restricting the critical frequency.

Stability of networks with non-zero dendritic rise and decay time constants and no time delays

The teardrop shaped zone and the time-delay effects on stability described in the previous section produce similar results to those seen in other studies (Marcus and Westervelt, 1989; Jirsa and Ding, 2004; Feng et al., 2006). However, our model also incorporates dendritic rise and decay time constants. In the next two sections we describe the stability of structural brain networks with non-zero dendritic time constants.

We first investigate brain networks with dendritic time constants and no propagation time delays. In this case τ = 0, α ≠ 0, and β ≠ 0. Hence L(ω) ≠ 1 and

From Equation (26) the stability boundary is given by

$\begin{array}{cc}{\lambda }_r=1-[1+2\gamma (1/\alpha +1/\beta )+{\gamma }^2/(\alpha \beta )]{\varpi }^2+{\gamma }^2/(\alpha \beta ){\varpi }^4,& (27)\end{array}$

and

$\begin{array}{cc}{\lambda }_i=-[2+\gamma (1/\alpha +1/\beta )]\varpi +[\gamma (1/\alpha +1/\beta )+2{\gamma }^2/(\alpha \beta )]{\varpi }^3.& (28)\end{array}$

where λ_r = Reλ, Imλ = λ_i, and ϖ is real. In this case λ − (ϖ) = 0 only has a finite number of solutions since (ϖ) is a polynomial of degree four. From Equation (28) the values of ϖ_c and ω_c are given by

$\begin{array}{cc}{\varpi }_c^2=\frac{2\alpha \beta +\gamma (\alpha +\beta )}{2{\gamma }^2+\gamma (\alpha +\beta )},& (29)\end{array}$

and

$\begin{array}{cc}{\omega }_c=\sqrt{\frac{2\alpha \beta \gamma +{\gamma }^2(\alpha +\beta )}{\alpha +\beta +2\gamma }},& (30)\end{array}$

respectively. An equation similar to Equation (30) was previously derived to describe gamma resonances produced by a similar mechanism (Robinson, 2005). If αβ < γ² then Equation (30) implies ω_c < γ.

Unlike the case for time delays in the previous section, the (ϖ) curve only crosses the real axis once for real ϖ > 0; an example of such a (ϖ) curve is shown in Figure 3A. However, the region enclosed by (ϖ) for −ϖ_c ≤ ϖ ≤ ϖ_c is again the stability zone and defined by Equation (25). The corresponding stability zone for the curve in Figure 3A is shown in Figure 3B.

FIGURE 3

Figure 3. Stability boundaries and stability zones for networks with non-zero dendritic rise and decay time constants and no time delays. The curve (ϖ) for real ϖ ≥ 0 (A) and corresponding stability zone and boundary (B) for networks with α = 60 s⁻¹, β = 240 s⁻¹ and γ = 100 s⁻¹ and τ = 0 s. The stability zone is shaded gray and the dark region within this zone is the unit disk. The dot-dashed line is the boundary for 1/α = 1/β = τ = 0. (C) Stability boundaries for β = 4α and γ/α = 100, 0.1, and 1.7. The values of γ/α are written next to their corresponding boundary. (D) Stability boundaries for γ/α = 1.7 and β/α = 1, 4 (dashed curve), and 10. The dot-dashed line in (C,D) is the boundary for 1/α = 1/β = τ = 0 and the dark gray region in (C,D) is the unit disk.

The stability zone in Figure 3B contains the unit disk and has a similar teardrop shape to the zone in Figure 2B with the arms of the parabola given by Equation (18) bent inward. In terms of stability, this implies non-zero 1/α and 1/β have similar effects to a propagation time delay. This is consistent with previous work on the corticothalamic model and highlights the low-pass filter effect of L(ω) (Robinson et al., 1997, 2001a,b; Rennie et al., 2002). The effective time delay resulting from γ/α and γ/β can be obtained by solving

$\begin{array}{cc}{e}^{-i{\varpi }_c\gamma \tau }=(1-i{\varpi }_c\gamma /\alpha )(1-i{\varpi }_c\gamma /\beta ),& (31)\end{array}$

for τ.

In the remainder of this section we explore the effect of varying γ, α, and β on the stability zone. We illustrate these effects by setting β/α to a positive constant. Firstly, for fixed β/α and large γ/α, Equations (27, 29) imply

$\begin{array}{cc}{\lambda }_r\approx 1+{(\gamma /\alpha )}^2{\varpi }^2({\varpi }^2-1),& (32)\end{array}$

and

$\begin{array}{cc}{\varpi }_c^2\approx \alpha /\gamma ,& (33)\end{array}$

respectively. From Equation (33), ω_c ≈ 0 for large γ/α and substituting Equation (33) into (32) shows that λ^c_r → −∞ as γ/α → ∞. These results are illustrated in Figures 3C, 4.

FIGURE 4

Figure 4. Change in stability boundary real axis intercept and critical frequency as a function of γ/α. (A) λ^c_r and (B) ω_c as a function of γ/α for γ τ = 0 and β/α = 1 (solid lines), 4 (dotted lines), and 10 (dot-dashed lines). In (B) ω_c is determined with γ = 100 s⁻¹. The dotted lines represent the nominal γ/α from Table 1.

The change in the stability zone for fixed β/α and varying γ/α is shown in Figure 3C. All the stability zones contain the unit circle, and for γ/α = 0.1 and 1.7 the stability zone is contained within the parabolic zone defined by Equation (18). However, the zone for γ/α = 100 has expanded so that its boundary intersects the parabolic boundary at λ_r ≈ −11. The ϖ_c in Figure 3C for γ/α = 100, 0.1, and 1.7 are 0.15, 4.0, and 1.08, respectively. The corresponding λ^c_r are −66 (intersection not seen), −19, and −4.9. For γ/α = 1.7, the nominal value from Table 1, the critical frequency is ω_c/2π ≈ 17 Hz.

The effect of changing β/α on the stability zone, while γ/α remains fixed, is shown in Figure 3D. The ϖ_c (λ^c_r) in Figure 3D for β/α = 1 and β/α = 10 are 0.77 (λ^c_r = −4.3) and 1.26 (λ^c_r = −6.2), respectively. Note that if β/α < 1, then the values of α and β can be swapped and the results in Figure 3 are reproduced. This is because exchanging γ/α and γ/β in Equation (4) has no effect on the dynamics and stability of a network. This implies the smallest stability zone with the minimum λ^c_r occurs when α = β. Note that experimental measurements of dendritic time constants in the brain give β/α ≈ 4–10 (Robinson et al., 2003a, 2004), the upper range in Table 1, and hence, a larger stability zone for brain activity.

Figure 3 suggests that as β/α increases the stability zone expands in a similar way to decreasing τ. However, as τ → 0 the stability zone converges to the parabolic zone (Equation 18), this is not the case for γ/β → ∞ and fixed γ/α. If γ/α is fixed to a positive constant and γ/β » γ/α then, from Equation (29), ϖ²_c ≲ 1 and therefore

$\begin{array}{cc}{\lambda }_r\approx 1-\gamma {\varpi }^2/\beta ,& (34)\end{array}$

and

$\begin{array}{cc}{\lambda }_i\approx -\gamma \varpi /\beta .& (35)\end{array}$

When λ_r = 0, Equation (34) implies $\varpi \approx \sqrt{\beta /\gamma }$ and ${\lambda }_i\approx -\sqrt{\gamma /\beta }$. Therefore as γ/β → ∞, the imaginary axis intercepts converge to ± ∞ and the stability zone expands to cover the entire region of the complex plane defined by Reλ < 1. Note that in this case, even though the eigenvalues can lie anywhere to the left of Reλ = 1, the dispersion solutions have an angular frequency ω < γ. These results explain the intersection of the stability boundary with the parabolic boundary in Figure 3C.

With Figure 3C these convergence results suggest that as γ/α increases from 0, the stability zone contracts toward the unit circle, and then expands again. Figure 4 shows the values of λ^c_r and ω_c as a function of γ/α and β/α. As γ/α increases from 0, λ^c_r rapidly increases from −∞ to a maximum value and then slowly decreases back to −∞; this decrease is greatest for the β/α = 1 curve which intersects the other two curves in Figure 4A. The maximum turning point for λ^c_r occurs when γ ≈ α for each β/α, with maximum λ^c_r decreasing as β/α increases. When β/α = 4 (which is the nominal value in Table 1), λ^c_r is approximately constant for 1 ≲ γ/α ≲ 3, with a maximum at γ/α ≈ 1.7, the nominal value. The corresponding ω_c curves in Figure 4B all show similar monotonic decreases as γ/α increases. The curves do not intersect or have a turning point as in Figure 4A. This decrease from ∞ is initially very rapid, before gradually decreasing to 0 as γ/α → ∞. This change approximately occurs at the nominal γ/α value in Table 1, where ω_c/2π ≲ 20 Hz for each β/α. Note that when β/α is a fixed constant, Equation (33) shows that ω_c → 0. Figure 4B also shows that increasing γ/α and decreasing γ/β results in a decreased value for ω_c.

In this section we have shown that physiologically realistic dendritic time constants have an effect on network stability similar to that of propagation time delays restricting, the critical frequency and the stability zone to a teardrop-shaped zone in the complex plane. However, unlike τ, for particular values of γ/α and γ/β the stability zone can expand to enclose an area outside the parabolic region described by Equation (18).

Stability of networks with time delays and non-zero dendritic time constants

The effect of having both time delays and non-zero dendritic time constants on stability is now described. In this case each of the parameters γ, 1/α, 1/β, and τ are non-zero and is given by Equation (11). As in the previous section, the stability boundary for these networks is defined by λ_r = Re(ϖ) and λ_i = Im(ϖ) for −∞ < ϖ < ∞. Analysis of these equations gives the properties of the stability boundary and zone. However, the effects on stability of having equal time delays and dendritic time constants are easily understood qualitatively as a combination of our previous results. Beginning with a network that has γ/α = γ/β = γτ = 0 and a parabolic stability zone, given by Equation (18), adding a time delay τ contracts the stability zone toward the unit circle, by “pulling in” the parabolic boundary, forming a teardrop-shaped zone within the original parabolic region. Adding dendritic rise and decay constants then, depending on their value, expands or contracts this stability zone. In all cases the stability zone is defined by

for 0 ≤ ϖ ≤ ϖ_c. Note that for large γτ, γ/α, and γ/β the contraction caused by γτ is greater than the expansion effects due to γ/α and γ/β because of the exponential in Equation (11).

In Figure 5 the stability zone of a brain network with plausible time delays and dendritic time constants (from Table 1) is shown. This shows that for realistic parameter values brain networks have a teardrop-shaped stability zone completely within the parabolic zone (Equation 18). For large realistic τ the stability zone is only slightly larger than the unit disc. The ϖ_c in Figure 5 for increasing γτ are 1.02, 0.72, and 0.35, respectively. The corresponding critical frequencies ω_c/2π for γ = 100 are then 16, 11.3, and 5.5 Hz, respectively.

FIGURE 5

Figure 5. Stability boundaries for brain networks with γ/α = 1.7, β/α = 4 and γτ = 0.1, 1.0, and 5.0. The values of γτ are written slightly to the left of their corresponding boundary. The gray region is the unit disk.

Stability and Dispersion Solutions of Randomly Connected Networks

We now investigate the stability and dispersion solutions of the randomly connected structural brain networks defined in the “Methods” section. For these networks we fix the model parameters γ, α, β, and τ to their nominal values in Table 1. The stability and dispersion solutions of this networks are determined from Sp(G). If all the eigenvalues λ in Sp(G) satisfy Equation (25) then the network is stable. The corresponding dispersion solutions are obtained numerically by solving Equation (10) for each λ using CROOT (Botten et al., 1983) (as described in “Methods”).

The spectrum of a RENs consists of one eigenvalue at npμ_e with the other n − 1 eigenvalues uniformly distributed in a disc of radius ${\mu }_e\sqrt{np(1-p)}<np{\mu }_e$. The spectrum of RCNs and RPNs with maximum μ_e and μ_i allowed by stability is distributed within the unit disk, see (Gray and Robinson, 2009a), with multiple eigenvalues near the stability boundary. Therefore, stability constrains the spectrum of random brain networks to the unit disk and the stability of random brain networks is independent of the γ, α, β, and τ. However, the frequencies of the dispersion solutions do depend on the model parameters.

In Figure 6 the spectrum and dispersion solutions for a REN, RCN, and RPN with the nominal model parameters are shown. The parameters of each network are set so that the networks are marginally stable and μ_e and |μ_i| are as large as possible while maintaining stability. Note the larger values of μ_e and |μ_i| for the RPN, compared to the RCN. This highlights that RPNs can have larger μ_e and |μ_i|, and hence be more responsive, before becoming almost certainly unstable as shown in (Gray and Robinson, 2009a). This suggests stability may have an effect on the arrangement of inhibitory and excitatory neurons and their physiology in structural brain networks.

FIGURE 6

Figure 6. Spectrum and dispersion solutions for a REN, a RCN, and a RPN. Each network has n = 50, p = 0.5, and the nominal parameters in Table 1, other parameters are set so that the network is close to marginal stability. The left column is the spectrum (crosses), the middle column shows the dispersion solutions (diamonds), and the right column shows an expanded view of the dispersion solutions near the imaginary axis with dotted lines representing the critical frequencies ± ω_c/2π. (A) REN with μ_e = 0.04, (B) RCN with p_i = 0.5, μ_e = −μ_i = 0.2, σ_e = σ_i = 0, and (C) RPN with p_i = 0.5, μ_e = −μ_i = 0.22, σ_e = σ_i = 0.

Each network has an infinite number of dispersion solutions because τ ≠ 0. The second column of Figure 6 shows the dispersion spectrum is symmetrically placed around the real axis. Each network has a qualitatively similar dispersion spectrum, with a finite cluster of solutions near the origin and a broad “arrowhead” of solutions for Imω ≲ 800 s⁻¹; this arrowhead has an infinite number of solutions with Imω decreasing as |Reω| increases. The third column shows the dispersion solutions near the imaginary axis. In Figure 6A3 one solution, ω = ω₁, lies on the origin, separate from the other dispersion solutions. This solution corresponds to the eigenvalue at λ ≈ 1 that is separated from the rest of the spectrum in Figure 6A1. This implies the dynamics of an REN will be dominated by a zero-frequency mode. Since the rest of the dispersion solutions have Imω « 0 all other modes rapidly decay to zero amplitude. However, the RCN and RPN in Figure 6 have a very similar dispersion spectrum with multiple dispersion solutions near the imaginary axis. The solutions closest to the imaginary axis have small frequencies ≲ 5 Hz. This shows that the presence of inhibitory connections allows random networks to have multiple marginally stable low frequency modes (Gray and Robinson, 2008, 2009a,b).

Discussion

We increased the physiological realism of a structural brain network model we studied previously in (Gray and Robinson, 2006, 2008, 2009a,b; Robinson et al., 2009) by allowing the network to have equal time delays τ for propagation between neuronal populations and non-zero dendritic rise 1/β and decay 1/α time constants. Under these assumptions the stability of an arbitrarily connected network of neural populations is determined by the network's gain matrix. The addition of time delays changed the stability zone in the complex plane from a parabolic region to a teardrop-shaped zone, dependent on α, β, τ, and the temporal damping rate γ. Our results are similar to previous work investigating the effect of time delays on stability of electrical activity within spatially continuous networks of neural tissue (Marcus and Westervelt, 1989; Jirsa and Ding, 2004; Feng et al., 2006; Qubbaj and Jirsa, 2007, 2009; Jirsa, 2009). This previous work has generally used integro-differential neural field equations with connectivity within a neural mass described by homogeneous or heterogeneous kernels. While in principle this work could be applied to large-scale connection topologies of discrete neural masses, as we have investigated here, this has generally not been done as it is difficult to incorporate arbitrary connectivity patterns (Qubbaj and Jirsa, 2009). In this work we have investigated the temporal dynamics of the overall electrical activity of arbitrarily connected structural brain networks, ignoring the spatial spread and propagation of electrical activity within individual neuronal populations.

In terms of stability the effect of non-zero dendritic time constants is similar to a time delay, further suggesting that dendrites act as a low-pass filter on synaptic inputs (Robinson et al., 1997, 2001a,b; Rennie et al., 2002). However, dendritic time constants can change the shape of the stability zone even allowing it to expand and enclose areas outside the parabolic stability region for networks with zero time delays and instantaneous dendritic rise and decay times. For all values of γ, α, β, and τ the stability zone contains the unit disk. This result implies that the stability criteria originally derived by May (1972, 1974) (that a network is stable if all its eigenvalues lie in the unit disk) is a sufficient condition for the stability of structural brain networks.

We also explored the dispersion solutions and frequencies of structural brain networks. If an initially stable brain network becomes unstable through a change in its connection gains, then at least one eigenvalue moves across the stability boundary and the network has an instability at a frequency given by the eigenvalue's corresponding dispersion solution. For networks with time delays and non-zero dendritic time constants there is a maximum frequency, the critical frequency ω_c/2π, at which initially stable networks will become unstable. For example, if all the gain matrix eigenvalues of a networks are initially inside the stability zone but then move across the stability boundary (e.g., due to changes in connection gains) then the frequency of these instabilities will be less than the critical frequency.

Measurements of brain activity (Stam et al., 1999; Robinson et al., 2001b; Breakspear, 2002; Breakspear et al., 2003) suggest the brain operates near marginal stability allowing the brain to have rich dynamics and a wide range of complex behavior. A network near marginal stability has eigenvalues near the stability boundary with corresponding modes that are the slowest to decay back to the steady state dominating the network's dynamics. These modes have a frequency less than the critical frequency. Using physiologically plausible parameter values in our structural network model (see Table 1) we would expect the electrical dynamics to be dominated by frequencies ≲100 Hz. When the nominal parameter values of α, β, γ, and τ are used the critical frequency is approximately 10 Hz and decreases as γτ increases (see Figures 4, 6).

For the randomly connected structural brain networks we investigated previously (Gray and Robinson, 2006, 2008, 2009a,b) the spectrum of the gain matrix is almost certainly contained in a disk centered on the origin with a radius dependent on the network's architecture and the average values of its excitatory and inhibitory gains. Thus, the stability zone of these networks is the unit disk and their stability is independent of time delays and dendritic time constants. Therefore, the results of that work remain valid in the more general case studied here. However, for the critical frequency is dependent on dendritic time constants, temporal damping rate, and time delays. We showed marginally stable randomly connected networks with inhibitory connections have multiple marginally stable low frequency dispersion solutions.

The primary goal of this and our previous work on structural brain networks is to understand how stability potential constrains the structure physiology of networks. For the randomly connected networks studied in (Gray and Robinson, 2006, 2008, 2009a,b), we have shown that time delays and non-zero dendritic time constants have minimal effect on their stability. One network type whose stability could be affected by these physiological properties is networks with inhibitory self-connections. The spectrum of these networks, even if they are randomly connected, is no longer restricted to a disc by stability but can have eigenvalues distributed within the teardrop-shaped region. Such networks could have marginally stable modes with frequencies (up to the critical frequency) in the alpha, beta, and gamma ranges. This will be explored in future work.

Likely Effects of Distributed or Varying Parameters on Stability

Assuming structural brain networks have equal γ, α, and β for each neural population and equal τ for each connection is unrealistic. Different neuronal populations in the brain have different parameter values; for example, excitatory cortical neurons have γ ≈ 100 s⁻¹ while for inhibitory cortical neurons γ ≈ 10³–10⁴ s⁻¹ (Robinson et al., 2004). Also, the time delay in real cortical networks is expected to vary from τ = 0 for self-connections to a large value for areas physically far apart. A realistic model of a structural brain network would therefore allow the model parameters γ, α, β, and τ to vary across neural populations. This variation could possibly be represented as a distribution.

The effect of distributed time delays on network stability has been studied using general models for network activity (Yi and Tan, 2002; Atay, 2003; Jirsa and Ding, 2004; Feng et al., 2006). In (Feng et al., 2006) and (Jirsa and Ding, 2004) networks with a distributed time delay with mean τ were shown to have a stability zone that contained the stability zone of networks with an constant delay equal to τ. These results are applicable to the brain networks studied here, since our model without dendritic time constants can be described as a specific case of the model studied in (Jirsa and Ding, 2004; Feng et al., 2006). This implies brain networks with distributed delay are more stable than networks with equal time delays; in the sense that a stable network with a distribution of delays could be unstable if its delays were replaced with a constant delay equal to the distribution mean. Hence, the equal time delay case is the least stable case and yields a bound on the stability of a structural brain network.

Here we have shown that plausible dendritic time constants have similar effects on stability as a time delay. This suggests that similar results to those found in (Jirsa and Ding, 2004; Feng et al., 2006) will likely be observed for distributed α and β. Also, we have shown that γ only affects the dispersion frequencies of a network, not its stability. A distributed γ is therefore expected to have no effect on the stability, only on the dispersion frequencies but this needs to be confirmed numerically.

Given the previous results investigating distributed parameter values we argue our results are still informative. But the exact effect of distributed γ, α, and β on structural brain network dynamics needs to be determined, particularly to understand the critical frequency and the dynamics of marginally stable modes in a networks electrical activity. To fully understand the stability and dynamics of structural brain networks with varying time delays and dendritic time constants requires a numerical approach. This will be investigated using CROOT (Botten et al., 1983) in future work.

Conclusions

We investigated the stability of discrete networks of neuronal populations using a simplified physiologically-based mean-field model of brain electrical activity. Incorporating time delays and non-zero dendritic time constants affects the stability of arbitrarily connected structural brain networks by constraining the eigenvalues of the gain matrix to a teardrop-shaped region in the complex plane. The stability of randomly connected networks of excitatory and inhibitory neuronal populations is unaffected by time delays and dendritic time constants; as stability constrains the gain matrix eigenvalues to the unit circle. However, the dispersion frequencies of instabilities are affected by network physiology. Randomly connected brain networks with the largest average excitatory and inhibitory gains allowed by stability can have multiple marginally stable low-frequency modes. Such networks would be highly responsive and adaptable to external stimuli while remaining stable, and have a wide range of flexible, adaptable, and complex behavior.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Acebrón, J. A., Bonilla, L. L., Vicente, C. J. P., Ritort, F., and Spigler, R. (2005). The kuramoto model: a simple paradigm for synchronization phenomena. Rev. Mod. Phys. 77, 137–185.

Allesina, S., and Tang, S. (2012). Stability criteria for complex ecosystems. Nature 483, 205–208.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Atay, F. M. (2003). Distributed delays facilitate amplitude death of coupled oscillators. Phys. Rev. Lett. 91:094101. doi: 10.1103/PhysRevLett.91.094101

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Atay, F. M., and Hutt, A. (2004). Stability and bifurcations in neural fields with finite propagation speed and general connectivity. SIAM J. Appl. Math. 65, 644–666.

Attwell, D., and Laughlin, S. B. (2001). An energy budget for signaling in the grey matter of the brain. J. Cereb. Blood Flow Met. 21, 1133–1145.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bassett, D. S., and Bullmore, E. (2006). Small-world brain networks. Neuroscientist 12, 1–12.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bechhoefer, J. (2005). Feedback for physicists: a tutorial essay on control. Rev. Mod. Phys. 77, 783–836.

Bollobás, B. (1985). Random Graphs. London: Academic Press.

Botten, L. C., Craig, M. S., and McPhedran, R. C. (1983). Complex zeros of analytic functions. Comput. Phys. Commun. 29, 245–259.

Breakspear, M. (2002). Nonlinear phase desynchronization in human electroencephalographic data. Hum. Brain Mapp. 15, 175–198.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Breakspear, M., Roberts, J. A., Terry, J. R., Rodrigues, S., Mahant, N., and Robinson, P. A. (2006). A unifying explanation of primary generalized seizures through nonlinear brain modeling and bifurcation analysis. Cereb. Cortex 16, 1296–1313.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Breakspear, M., Terry, J. R., and Friston, K. J. (2003). Modulation of excitatory synaptic coupling facilitates synchronization and complex dynamics in a nonlinear model or neuronal dynamics. Neurocomputing 52–54, 151–158.

Bressloff, P. (2011). Spatiotemporal dynamics of continuum neural fields. J. Phys. A Math. Theor. 45:033001. doi: 10.1088/1751-8113/45/3/033001

CrossRef Full Text

Bullmore, E., and Sporns, O. (2009). Complex brain networks: graph theoretical analysis of structural and functional systems. Nat. Rev. Neurosci. 10, 186–198.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Bullmore, E., and Sporns, O. (2012). The economy of brain network organization. Nat. Rev. Neurosci. 13, 336–349.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Coombes, S. (2005). Waves, bumps, and patterns in neural field theories. Biol. Cybern. 93, 91–108.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Coombes, S., Venkov, N. A., Shiau, L., Bojak, I., Liley, D. T., and Laing, C. R. (2007). Modeling electrocortical activity through improved local approximations of integral neural field equations. Phys. Rev. E 76:051901. doi: 10.1103/PhysRevE.76.051901

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Cvetković, D., Doob, M., and Scahs, H. (1995). Spectra of Graphs. 3rd Edn. Heidelberg: Johann Ambrosius Barth Verlag.

da Silva, F. H. L., Hoeks, A., Smits, A., and Zetterberg, L. H. (1974). Model of brain rhythmic activity: the alpha-rhythm of the thalamus. Kybernetik 15, 27–37.

Pubmed Abstract | Pubmed Full Text

Deco, G., Jirsa, V., Robinson, P., Breakspear, M., and Friston, K. (2008). The dynamic brain: from spiking neurons to neural masses and cortical fields. PLoS Comput. Biol. 4:e1000092. doi: 10.1371/journal.pcbi.1000092

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Felleman, D. J., and van Essen, D. C. (1991). Distributed hierarchical processing in the primate cerebral cortex. Cereb. Cortex 1, 1–47.

Pubmed Abstract | Pubmed Full Text

Feng, J., Jirsa, V. K., and Ding, M. (2006). Synchronization in networks with random interactions: theory and applications. Chaos 16:015109. doi: 10.1063/1.2180690

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Freeman, W. J. (1975). Mass Action in the Nervous System: Examination of the Neurophysiological Basis of Adaptive Behavior Through the EEG. New York, NY: Academic Press.

Gray, R., and Robinson, P. (2006). Stability and spectra of randomly connected excitatory cortical networks. Neurocomputing 70, 1000–1012.

Gray, R., and Robinson, P. (2008). Stability and synchronization of random brain networks with a distribution of connection strengths. Neurocomputing 71, 1373–1387.

Gray, R. T., and Robinson, P. A. (2009a). Stability and structural constraints of random brain networks with excitatory and inhibitory neural populations. J. Comput. Neurosci. 27, 81–101.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Gray, R. T., and Robinson, P. A. (2009b). Stability of random brain networks with excitatory and inhibitory connections. Neurocomputing 72, 1849–1858.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hilgetag, C.-C., Burns, G. A. P. C., O'Neill, M. A., Scannell, J. W., and Young, M. P. (2000a). Anatomical connectivity defines the organization of clusters of cortical areas in the macaque monkey and the cat. Philos. Trans. R. Soc. Lond. B 355, 91–110.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hilgetag, C.-C., O'Neill, M. A., and Young, M. P. (2000b). Hierarchical organization of macaque and cat cortical sensory systems explored with a novel network processor. Philos. Trans. R. Soc. Lond. B 355, 71–89.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Hogg, T., Huberman, B. A., and McGlade, J. M. (1989). The stability of ecosystems. Proc. R. Soc. B 237, 43–51.

Horn, R. A., and Johnson, C. R. (1985). Matrix Analysis. New York, NY: Cambridge University Press.

Jirsa, V. (2009). Neural field dynamics with local and global connectivity and time delay. Philos. Trans. R. Soc. Lond. A Math. Phys. Eng. Sci. 367, 1131–1143.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jirsa, V. K., and Ding, M. (2004). Will a large complex system with time delays be stable? Phys. Rev. Lett. 93:070602, 1–4. doi: 10.1103/PhysRevLett.93.070602

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jirsa, V. K., and Haken, H. (1996). Field theory of electromagnetic brain activity. Phys. Rev. Lett. 77, 960–963.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Jouve, B., Rosenstiehl, P., and Imbert, M. (1998). A mathematical apporach to the connectivity between the cortical visual areas of the macaque monkey. Cereb. Cortex 8, 28–39.

Pubmed Abstract | Pubmed Full Text

Kandel, E. R., Schwartz, J. H., and Jessell, T. M. (2000). Principles of Neural Science. 4th Edn. New York, NY: McGraw-Hill.

Kim, J. W., and Robinson, P. A. (2007). Compact dynamical model of brain activity. Phys. Rev. E 75:031907. doi: 10.1103/PhysRevE.75.031907

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Koch, C. (2004). The Quest for Consciousness: A Neurobiological Approach. Englewood: Roberts and Company.

Kuramoto, Y. (1984). Chemical Oscillations, Waves, and Turbulence. New York, NY: Dover.

Laughlin, S. B., de Ruyter van Stevenick, R. R., and Anderson, J. C. (1998). The metabolic cost of neural information. Nat. Neurosci. 1, 36–41.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Lennie, P. (2003). The cost of cortical computation. Cogn. Brain Res. 13, 493–497.

Pubmed Abstract | Pubmed Full Text

Marcus, C. M., and Westervelt, R. M. (1989). Stability of analog neural networks with delay. Phys. Rev. A 39, 347–359.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

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

Pubmed Abstract | Pubmed Full Text

May, R. M. (1974). Stability and Complexity in Model Ecosystems. 2nd Edn. Princeton, NJ: Princeton University Press.

McCann, K. S. (2000). The diversity-stability debate. Nature 405, 228–233.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Mitrinović, D. S., and Keckić, J. D. (1984). The Cauchy Theory of Residues: Theory and Applications. Dordrecht: D. Reidel Publishing Company.

Murray, J. D. (2002). Mathematical Biology. 3rd Edn. New York, NY: Springer.

Nunez, P. L. (1974). Wave-like properties of the alpha rhythm. IEEE Trans. Biomed. Eng. 21, 473–482.

Nunez, P. L. (1995). Neocortical Dynamics and Human EEG Rhythms. New York, NY: Oxford University Press.

Pecora, L. M., and Carroll, T. L. (1998). Master stability functions for synchronized coupled systems. Phys. Rev. Lett. 80, 2109–2112.

Qubbaj, M., and Jirsa, V. (2007). Neural field dynamics with heterogeneous connection topology. Phys. Rev. Lett. 98:238102. doi: 10.1103/PhysRevLett.98.238102

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Qubbaj, M., and Jirsa, V. (2009). Neural field dynamics under variation of local and global connectivity and finite transmission speed. Physica D 238, 2331–2346.

Rennie, C. J., Robinson, P. A., and Wright, J. J. (2002). Unified neurophysical model of EEG spectra and evoked potentials. Biol. Cybern. 86, 457–471.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Robinson, P. A. (2003). Neurophysical theory of coherence and correlations of electroencephalographic and electrocorticographic signals. J. Theor. Biol. 222, 163–175.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Robinson, P. A. (2005). Propagator theory of brain dynamics. Phys. Rev. E 72:0119404. doi: 10.1103/PhysRevE.72.011904

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Robinson, P. A., Henderson, J. A., Matar, E., Riley, P., and Gray, R. T. (2009). Dynamical reconnection and stability constraints on cortical network architecture. Phys. Rev. Lett. 103, 1–4.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Robinson, P. A., Loxley, P. N., O'Connor, S. C., and Rennie, C. J. (2001a). Modal analysis of corticothalamic dynamics, electroencephalographic spectra, and evoked potentials. Phys. Rev. E 63:041909. doi: 10.1103/PhysRevE.63.041909

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Robinson, P. A., Rennie, C. J., Wright, J. J., Bahramali, H., Gordon, E., and Rowe, D. L. (2001b). Prediction of electroencephalographic spectra from neurophysiology. Phys. Rev. E 63:021903. doi: 10.1103/PhysRevE.63.021903

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Robinson, P. A., Rennie, C. J., and Rowe, D. L. (2002). Dynamics of large-scale brain activity in normal arousal states and epileptic seizures. Phys. Rev. E 65:041924. doi: 10.1103/PhysRevE.65.041924

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Robinson, P. A., Rennie, C. J., Rowe, D. L., and O'Connor, S. C. (2004). Estimation of multiscale neurophysicologic parameters by electroencephalographic means. Hum. Brain Mapp. 23, 53–72.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Robinson, P. A., Rennie, C. J., Rowe, D. L., O'Connor, S. C., Wright, J. J., Gordon, E., et al. (2003a). Neurophysical modeling of brain dynamics. Neuropsychopharmacology 28, S74–S79.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Robinson, P. A., Whitehouse, R. W., and Rennie, C. J. (2003b). Nonuniform corticothalamic continuum model of electroencephalographic spectra with application to split-alpha peaks. Phys. Rev. E 68:021922. doi: 10.1103/PhysRevE.68.021922

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Robinson, P. A., Rennie, C. J., and Wright, J. J. (1997). Propagation and stability of waves of electrical activity in the cerebral cortex. Phys. Rev. E 56, 826–840.

Robinson, P. A., Rennie, C. J., Wright, J. J., and Bourke, P. D. (1998). Steady states and global dynamics of electrical activity in the cerebral cortex. Phys. Rev. E 58, 3557–3571.

Rubinov, M., and Sporns, O. (2010). Complex network measures of brain connectivity: uses and interpretations. Neuroimage 52, 1059–1069.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Scannell, J. W., Blakemore, C., and Young, M. P. (1995). Analysis of connectivity in the cat cerebral cortex. Int. J. Neural Syst. 15, 1463–1483.

Pubmed Abstract | Pubmed Full Text

Sporns, O. (2011). The human connectome: a complex network. Ann. N.Y. Acad. Sci. 1224, 109–125.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sporns, O., Chialvo, D. R., Kaiser, M., and Hilgetag, C. C. (2004). Organization, development and function of complex brain networks. TRENDS Neurosci. 8, 418–425.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sporns, O., Tononi, G., and Edelman, G. M. (2000). Theoretical neuroanatomy: relating anatomical and functional connectivity in graphs and cortical connection matrices. Cereb. Cortex 10, 127–141.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sporns, O., Tononi, G., and Kötter, R. (2005). The human connectome: a structural description of the human brain. PLoS Comput. Biol. 1:e42. doi: 10.1371/journal.pcbi.0010042

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Sporns, O., and Zwi, J. D. (2004). The small world of the cerebral cortex. Neuroinformatics 2, 145–161.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Stam, C., Pijn, J. P. M., Suffczynski, P., and da Silva, F. H. L. (1999). Dynamics of the human alpha rhythm: evidence for non-linearity? Clin. Neurophysiol. 110, 1801–1813.

Pubmed Abstract | Pubmed Full Text

Steriade, M., Gloor, P., Linas, R. P., da Silva, F. H. L., and Mesulam, M. M. (1990). Basic mechanisms of cerebral rhythmic activities. Electroencephalogr. Clin. Neurophysiol. 76, 481–508.

Pubmed Abstract | Pubmed Full Text

Steuer, R. (2007). Computational approaches to the topology, stability and dynamics of metabolic networks. Phytochemistry 68, 2139–2151.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Taverna, D. M., and Goldstein, R. A. (2002). Why are proteins marginally stable? Proteins 46, 105–109.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Turing, A. M. (1952). The chemical basis of morphogenesis. Philos. Trans. R. Soc. Lond. B 237, 37–72.

van den Heuvel, M., and Sporns, O. (2011). Rich-club organization of the human connectome. J. Neurosci. 31, 15775–15786.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Venkov, N. A., Coombes, S., and Matthews, P. C. (2007). Dynamic instabilities in scalar neural field equations with space-dependent delays. Physica D 232, 1–15.

Wen, Q., and Chklovskii, D. B. (2005). Segregation of the brain into gray and white matter: a design minimizing conduction delays. PLoS Comput. Biol. 1:e78. doi: 10.1371/journal.pcbi.0010078

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Wilson, H. R., and Cowan, J. D. (1973). A mathematical theory of the functional dynamics of cortical and thalamic nervous tissue. Kybernetik 13, 55–80.

Pubmed Abstract | Pubmed Full Text

Wright, J. J., and Liley, D. T. J. (1996). Dynamics of the brain at global and microscopic scales: neural networks and the EEG. Behav. Brain Sci. 19, 285–320.

Wright, J. J., Robinson, P. A., Rennie, C. J., Gordon, E., Bourke, P. D., Chapman, C. L., et al. (2001). Towards an integrated continuum model of cerebral dynamics: the cerebral rhythms, synchronous oscillation and cortical stability. Biosystems 63, 71–88.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Yi, Z., and Tan, K. K. (2002). Dynamic stability conditions for Lotka-Volterra recurrent neural networks with delay. Phys. Rev. E 66:011910, 1–8. doi: 10.1103/PhysRevE.66.011910

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text

Young, M. P. (2000). The architecture of visual cortex and inferential processes in vision. Spat. Vis. 13, 137–146.

Pubmed Abstract | Pubmed Full Text | CrossRef Full Text