## Abstract

Spatially continuous networks with heterogeneous connections are ubiquitous in biological systems, in particular neural systems. To understand the mutual effects of locally homogeneous and globally heterogeneous connectivity, we investigate the stability of the rest-state activity of a neural field as a function of its connectivity. The variation of the connectivity is operationalized through manipulation of a heterogeneous two-point connection embedded into the otherwise homogeneous connectivity matrix, as well as by variation of connectivity strength and a finite transmission speed. The latter results in a time delay of communication among individual brain areas. We demonstrate that the local connectivity generates the well-known power-law behaviour of the electroencephalographic power spectrum with an exponent close to −2, whereas the global connections generate a more characteristic line spectrum. These spectral characteristics are routinely observed in large-scale topographies of the human brain.

## 1. Introduction

Large-scale neural network models are thought to be involved in the implementation of cognitive function of the brain (Bressler 1995, 2002, 2003; Bullmore *et al*. 1996; Mesulam 1998; Mountcastle 1998; McIntosh 2000; Bressler & Kelso 2001; Jirsa 2004; Bressler & Tognoli 2006; Bressler & McIntosh 2007). In particular, peripheral areas are more functionally differentiated, whereas the ensuing cognitive integration appears to require more global network activations. The properties of the global network dynamics will then naturally depend on the large-scale, i.e. global, connectivity of the network components, as well as their local connectivity and dynamics (Sporns & Tononi 2002, 2007; Sporns 2003; Jirsa 2004; Beggs *et al*. 2007). So far, theoretical efforts have focused almost exclusively on the study of networks with discretely connected nodes and complicated connectivity, but simple local dynamics, in which all neural activity is lumped into a single neural mass (Lopes da Silva *et al*. 1974; Freeman 1975; van Rotterdam *et al*. 1982; Tagamets & Horwitz 1998; David & Friston 2003; Jirsa & Ding 2004; Campbell *et al*. 2006; see Ermentrout (1998) and Jirsa (2004) for reviews).

As an additional complication, large-scale brain networks will be affected by time delays via signal transmission along the connecting pathways. The time delays may reach up to 200 ms (Nunez 1995) for the human brain, which is on a similar time scale as human brain function and hence is not negligible. The space–time structure of the couplings, i.e. coupling strength between areas (space) and the time delay via transmission between areas (time), has specifically been shown to contribute to the emergence of spontaneous coherent fluctuations in the resting brain (Ghosh *et al*. 2008*a*,*b*). Here, the authors used a biologically realistic primate connectivity matrix with time delays. The individual network nodes were neural population oscillators embedded into a network resting in its equilibrium state. Under the influence of noise, the network dynamics intermittently exhibited damped wave propagation, which was associated with the spatiotemporal waves observed in human brain imaging studies. The authors demonstrated that spatiotemporal correlations observed in functional magnetic resonance imaging are captured by the network model only when time delays are taken into account. These correlations reflect correlated oscillatory activity on ultra-slow scales, less than 0.1 Hz. On faster time scales, 10–500 ms, the neural activity displays irregular, though characteristic, rhythms captured by electroencephalography (EEG) and magnetoencephalography (MEG). Although the modelling study by Ghosh *et al*. (2008*a*,*b*) explains the generation of these rhythms and the dominant peaks in the power spectra, it entirely lacks the ubiquitous characteristic power-law behaviour of spectral power (see Buzsaki (2006) for a range of exponents). A possible explanation for the absence of this feature is the implementation of single neural oscillators as the network nodes rather than the more realistic spatially expanded neural fields of population activity. Neural fields describe the temporal change of neural activity on a local spatial scale, typically within a brain area (Wilson & Cowan 1972; Wilson *et al*. 1973; Amari 1977; Bressloff & Coombes 1997; Atay & Hutt 2005; Hutt & Atay 2005; see Ermentrout (1998), Jirsa (2004) and Coombes (2005) for reviews). These approaches use a translationally invariant, so-called homogeneous, connectivity and also take time delays into consideration, which, however, given the small spatial extent of a brain area (as compared to the whole brain), play a smaller role.

When considering large-scale brain dynamics as in the resting brain activity, or generally for all cognitive function, theoretical and analytical means must be developed that are specifically targeted to the properties of large-scale network dynamics integrating local and global contributions. Existing attempts so far include neural field theories that approximate the large-scale components of the connectivity matrix as translationally invariant and exponentially decaying over space, but on a larger spatial scale than biologically realistic (Nunez 1974; Wright & Liley 1995; Jirsa & Haken 1996; Jirsa *et al*. 1997; Robinson *et al*. 1997; Breakspear 2004). These approaches have been successful in capturing various phenomena of large-scale brain dynamics, including characteristic EEG power spectra (Nunez 1974; Robinson *et al*. 1997), epilepsy (Breakspear *et al*. 2006) and MEG activity during sensorimotor coordination (Jirsa & Haken 1996; Jirsa *et al*. 1997). No systematic investigation, however, has been performed so far to test the validity of the translationally invariant approximation. There are only a few recent studies exploring the effects of global connecting pathways upon the local dynamics of coupled brain areas (see Jirsa & Kelso 2000; Qubbaj & Jirsa 2007). In the current paper, we wish to focus upon the contributions of the local connectivity to the large-scale dynamics of the entire network, in particular on its effects captured by the power spectrum. To do so, we identify a convenient mathematical representation for the network dynamics with a locally invariant (homogeneous) and globally variant (heterogeneous) architecture via a one-dimensional spatiotemporally continuous integro-differential field equation with space-dependent delay (following Jirsa & Kelso 2000; Qubbaj & Jirsa 2007). We demonstrate that the local connections contribute to the generation of spectral power-law behaviour with an exponent of −2, whereas the global connectivity is involved in the generation of the characteristic frequency peaks.

## 2. Mathematical basis of neural field dynamics

Let *ψ*(*x*,*t*) be the neural field capturing the population activity at time point *t* and position *x*. The dynamics of the neural field can then be described by the following integro-differential equation:(2.1)where 1/*ϵ*>0∈ represents the intrinsic time scale. The dot indicates the first time derivative. The spatial domain of the neural field is denoted by *Γ*, where *x*∈*Γ*=[0,*L*] and *L* is the spatial length of the neural field. The homogeneous connectivity function, *W*_{hom}(|*x*−*y*|), is translationally invariant. If the connectivity function does not have such properties, then we call it heterogeneous, *W*_{het}(*x*,*y*)≠*W*_{het}(|*x*−*y*|). The parameters *c* and *v* are the propagation velocities through the homogeneous and heterogeneous connections, respectively.

### (a) Local homogeneous connectivity is translationally invariant

For the moment, consider only the homogeneous connectivity term. By using the properties of delta functions, the system can be written as follows:(2.2)By using the Green's function method (Nunez 1995; Jirsa & Haken 1996; Jirsa *et al*. 1997), we write (2.2) as(2.3)where *G*(*x*−*y*,*t*−*T*)=*W*_{hom}(|*x*−*y*|)*δ*(*t*−*T*+|*x*−*y*|/*c*). Using Fourier transformations, *G*(*x*−*y*,*t*−*T*) can be written as(2.4)and it follows for equation (2.3) that(2.5)where *Ψ*(*k*,*ω*), and are the Fourier transformations of the functions *ψ*(*x*,*t*), *G*(*x*,*t*), and *S*(*x*,*t*), respectively. Apply the time derivative and rearrange the terms to obtain(2.6)The last two integrals on the right-hand side (r.h.s.) are equivalent to and , respectively. Hence, equation (2.6) reduces to(2.7)or simply(2.8)

For concreteness, we compute for the often used exponentially decaying connectivity function (Nunez 1974; Jirsa & Haken 1996; Jirsa *et al*. 1997):(2.9)where *σ* represents the spatial width of the distribution. Then we obtain(2.10)where *ξ*=*x*−*y* and *τ*=*t*−*T*.

For short-range interactions, *σ*≪*L*, the dispersion relation can be approximated as(2.11)giving rise to a diffusion term *σ*^{2}*k*^{2} for the lower orders of the expansion. Higher orders are indicated by the dots. The nearest-neighbour diffusive interaction is more evident when we return to the space–time domain and rewrite (2.1) as(2.12)where we absorbed the expression (1+*σ*/*c*) in the constants and *W*_{het}. If the homogeneous interactions are of longer range, then higher orders of the spatial derivatives will be involved, resulting in terms as present in the Swift–Hohenberg equation (including a fourth-order spatial derivative; Swift & Hohenberg 1977). The latter is known to give rise to self-sustained spatial patterns that are thought to be involved in working memory.

### (b) Neural fields with locally homogeneous and heterogeneous connectivity

Assume the heterogeneous two-point connection, *W*_{het}(*x*,*y*), as follows:(2.13)where *μ*_{ij}∈ represents the coupling strength through the heterogeneous connection. Now, we can write equation (2.1) as(2.14)Our main interest is the study of the stability of the rest state of a cortical architecture with local and global connectivity. Hence we consider the linear stability of the fixed point solution *ψ*_{0}(*x*) with . We perform a mode expansion of *ψ*(*x*,*t*) into a set of spatial basis functions {*ϕ*_{k}(*x*)},(2.15)and is the adjoint set satisfying the biorthogonality condition(2.16)where *δ*_{kl} is the Kronecker delta function. In (2.15) *ξ*_{k}(*t*) is the time-dependent amplitude related to the spatial basis function *ϕ*_{k}(*x*). Linearizing around *ψ*_{0}(*x*), and using (2.15) in (2.14), we obtain(2.17)Multiplying both sides by and integrating over the whole space, we obtain(2.18)with *d*=|*x*_{2}−*x*_{1}|. If we use the results given in (2.12), we can rewrite (2.14) as(2.19)Equation (2.19) describes the linearized dynamics of the neural field's equilibrium state in the presence of local homogeneous and global heterogeneous pathways, and for sufficiently short local path length *σ*. An illustration of its architecture is shown in figure 1.

## 3. Neural fields with local diffusion and global two-point heterogeneous connections

Consider a diffusively coupled neural field with a heterogeneous connection at two points of the field. Then the neural field dynamics is determined by equation (2.19). If we set *D*=0, the system reduces to that of two coupled oscillators described by the following delay differential equation:(3.1)with *i*,*j*=1,2, *i*≠*j*. The time-dependent amplitude of the *i*th oscillator is defined by *x*_{i}∈ and the discrete time delay by . Letting * B*=(

*μ*

_{ij}) be the coupling matrix, the system can be written in the vector/matrix form as(3.2)Following Jirsa & Ding (2004), we can decompose

**according to**

*B**=*

**B**

**E**^{−1}

*, where*

**ΛE***is the Jordan form with eigenvalues and matrix*

**Λ***contains the corresponding eigenvectors . Multiplying (3.2) from the left with*

**E***we obtain(3.3)If we set the eigenmodes to , equation (3.3) reduces to a decoupled representation of the dynamics of the system in terms of the eigenmodes as follows:(3.4)The stability properties of the above one-dimensional delay differential equation have been widely studied. Assuming a solution of the form*

**E***u*(

*t*)=e

^{λt},

*λ*∈, the stability condition is determined by the following characteristic equation:(3.5)Notice that any solution of (3.4) satisfying Re[

*λ*]<0 is stable and when Re[

*λ*] becomes positive the system destabilizes. On the critical surface parametrized by

*ϵ*,

*μ*

_{12},

*μ*

_{21}and

*τ*, the solution

*λ*of the characteristic equation is purely imaginary, i.e.,

*λ*=Im[

*λ*]=i

*ω*, . This surface represents the boundary at which the system changes stability and is shown in figure 2. A stability change through Re[

*λ*]=∞ is impossible (Datko 1978), hence every instability must occur at the surface defined by

*λ*=i

*ω*.

Using the characteristic polynomial *H*(*λ*) for *λ*=i*ω* with some algebra we can show that(3.6)where *K*=*βτ*+e^{iωτ}. The above result implies that, as *τ* increases across the critical surface from below, the equilibrium state always destabilizes and remains unstable for all larger *τ*.

The example of two coupled oscillators, *D*=0, represents the extreme case where only global interactions are present. To study the interplay between the local and global modes in the system with local diffusive coupling, we need to consider *D*≠0. Note that *D* must be sufficiently large to have an effect. In turn, this contributes significantly to the large-scale dynamics of the system, for which we assume the following ansatz:(3.7)Linearize (2.19) and insert (3.7) to obtain(3.8)Multiplying both sides by e^{−ikx} and integrating over space, we obtain(3.9)where we have made use of the fact that . To study the stability of (3.9), let us consider the simplest example where one mode will account for most of the dynamics of the system, i.e., *ψ*(*x*,*t*)=*ξk*(*t*)e^{ikx}+c.c., where c.c. denotes the complex conjugate. Hence, (3.9) reduces to(3.10)where *P*=*P*_{1}+i*P*_{2}, , and *d*=|*x*_{2}−*x*_{1}|=*x*_{2}−*x*_{1}, *x*_{2}>*x*_{1}. The characteristic equation of (3.10) is(3.11)Equation (3.11) characterizes the dynamics and provides the conditions for an instability of the system due to a heterogeneous connection with delay. Jirsa and colleagues (Jirsa & Ding 2004; Qubbaj & Jirsa 2007) discussed a similar characteristic equation with *ϵ*+*Dk*^{2}=1 and obtained the critical stability surface for the parameters *P*_{1}, *P*_{2} and *τ*. To get an insight into the essential nature of the characteristic equation (3.11) and the effect of diffusion on the stability surface, we will discuss the fully symmetric case, *μ*_{12}=*μ*_{21}=*μ*. The approach generalizes for any combination of *μ*_{12} and *μ*_{21} (see also Jirsa & Ding 2004; Qubbaj & Jirsa 2007). Then the characteristic equation (3.11) reduces to(3.12)Separating real and imaginary parts, with *λ*=i*ω* as earlier, we obtain(3.13)(3.14)Both equations must be satisfied for an instability to exist. Squaring and summing the two right-hand sides of (3.13) and (3.14), we obtain the dispersion relation *ω*^{2}=*μ*^{2}cos^{2}(*kd*)−(*Dk*^{2}+*ϵ*)^{2}, which is plotted in figure 3*a*. For increasing wavenumber *k*, fewer frequencies *ω* are available, since *ω*^{2} becomes negative. The smaller the diffusion constant *D* is, the more frequencies will be available. Consider equation (3.13), which is plotted in figure 3*b* for the critical value *ωτ*=0 (at which |cos(*ωτ*)|≤1 is the largest). There will be an intersection of the two curves that gives the value of *k* of the destabilizing mode. It is evident from figure 3*b* that the uniform mode, *k*=0, is the mode closest to the instability for *μ*>0 with *D*>0 and *ϵ*>0. If there are multiple intersections (depending on the time delay *τ* and other parameters), a discrete line spectrum will be possible. Modes with increasing wavenumber have an increasingly smaller degree of instability, when the vertical distance in figure 3*b* between *μτ* cos(*kd*) and *Dk*^{2}+*ϵ* increases. Since the degree of stability scales monotonically with the distance between these curves, our reasoning holds. In other words, modes with larger wavenumber *k* are more strongly damped.

We illustrate these effects in a simulation of equation (2.12) with white noise. The parameters are chosen such that the equilibrium state of the neural field is stable, but close to instability. As a consequence, the fluctuations continually kick the system out of its equilibrium, which then exhibits a characteristic transient dynamics back to equilibrium. The global field power, i.e., the spatial integral over the power spectra of the neural field, is plotted in figure 4 for various configurations of connectivity. In all cases there is a 1/*f*^{2} behaviour of the spectra for higher frequencies. However, the neural fields with homogeneous connectivity show enhanced power in the low-frequency regions, which is not the case for the purely heterogeneously connected network. The discrete spectral peak in figure 4 is a consequence of the interactions between the two heterogeneously coupled regions. To allow for a comparison with experimental data, we show the power spectrum of human rest-state EEG in figure 5. The data have been obtained from a single electrode located in the centre of the scalp (Cz) during a 4 min recording session with eyes closed. The power spectrum shows power-law behaviour with an exponent of slightly less than −2 over a large frequency range with a spectral peak in the so-called alpha range, approximately 8–12 Hz.

## 4. Conclusions

In this paper, we discussed a key issue in the field of biologically realistic large-scale networks, that is, the interplay between local and global architectures leading to a large-scale network dynamics. To gain insight into the interplay between local and global architectural components, we operationalized the global connectivity by a two-point connection with a finite transmission speed. In real networks, obviously more complicated connectivity prevails, resulting in a rich spatiotemporal brain dynamics on multiple scales (see Ghosh *et al*. (2008*b*) for a detailed modelling study with realistic primate connectivity). Evidently the resulting power spectra will also be a function of the degree of complexity of the large-scale connectivity skeleton. For this reason, our simplistic architecture and neural field dynamics may serve only as a toy model to discuss the interplay between local homogeneous and global heterogeneous brain connectivity effects. As a consequence, the details of the dynamic effects of certain parameter changes should not be thought of as representative for the fully connected brain.

For instance, our toy model of a two-point connection in a neural field shows oscillatory instabilities for *μ*_{12}=−*μ*_{21} (see figure 4), which corresponds to a link with an excitatory terminal on one end and an inhibitory terminal on the other. Such a configuration is not biologically realized, since all large-scale pathways have excitatory terminals. For purely excitatory fibres, *μ*_{12}, *μ*_{21}>0, the first instability is mostly non-oscillatory, as we showed in the previous section and in figure 3. Higher instabilities include oscillatory dynamics of non-uniform spatial modes, but occur when the spatially uniform mode, *k*=0, is already unstable. Still they are visible in the presence of noise, when the neural field returns to its stable rest state, but less prevalent. Although, it is evident that an arbitrary connectivity matrix can be constructed from the linear sum of two-point connections, an understanding of how the connectivity influences the network dynamics remains to be developed. It is in this sense that our discussion of a two-point connection and network dynamics should be seen as illustrative. However, it does shed light on the following aspect of large-scale brain dynamics: a heterogeneous large-scale connectivity provides the skeleton for a set of neural populations connected via pathways with transmission delays. It can be as simple as presented here, i.e. a two-point connection, or as biologically realistic (and complicated) as in Ghosh *et al*. (2008*b*). When the neural populations at the nodes of the network are expanded from a point in space to a neural field with local short-range connectivity in space, then this expansion in space will enhance the 1/*f*^{2} behaviour in the power spectrum, as we showed here. Since this characteristic is in the first place independent of the details of the local connectivity (as long as the rest state remains stable), the power-law behaviour with an exponent close to −2 will be ubiquitous.

## Acknowledgements

This research was funded by the grants Brain NRG JSM22002082, ATIP (CNRS) and the CNRS programme Neuroinformatics. We wish to thank Murad Qubbaj for preparing figure 2.

## Footnotes

- © 2009 The Royal Society