## Abstract

The study of dynamical systems defined on complex networks provides a natural framework with which to investigate myriad features of neural dynamics and has been widely undertaken. Typically, however, networks employed in theoretical studies bear little relation to the spatial embedding or connectivity of the neural networks that they attempt to replicate. Here, we employ detailed neuroimaging data to define a network whose spatial embedding represents accurately the folded structure of the cortical surface of a rat brain and investigate the propagation of activity over this network under simple spreading and connectivity rules. By comparison with standard network models with the same coarse statistics, we show that the cortical geometry influences profoundly the speed of propagation of activation through the network. Our conclusions are of high relevance to the theoretical modelling of epileptic seizure events and indicate that such studies which omit physiological network structure risk simplifying the dynamics in a potentially significant way.

## 1. Introduction

The newly emerging discipline of network science provides a general framework for representing, modelling and predicting the behaviour of complex systems belonging to areas as diverse as social science, biology and information technology [1]. Motivated by the observation that most real-world networks fail to conform to the homogeneous Poissonian degree distribution admitted by Erdo˝s & Rényi [2] random graphs, improved network models were constructed (most notably the *small-world model* of Watts & Strogatz [3] and the *preferential attachment model* of Barabási & Albert [4]) that were capable of recovering many of the interesting features displayed by real-world network data. Initial investigations into complex networks focused primarily on the characterization of networks in terms of a small number of topological parameters; however, more recently, interest has shifted towards understanding the influence of network structure on the dynamic processes occurring upon them—see [5,6] and references therein.

Such investigations are particularly relevant to biological systems, in which a well-defined network structure is frequently a key feature, the evolution and topology of which are presumed to affect relevant biological processes. A paradigm for such studies is that of neural systems, which have been widely studied in this context with the aim of providing a more complete understanding of epilepsy and other neural conditions [7]. Epilepsy is characterized by recurrent and unpredictable instances of ‘excessive or synchronous neuronal activity’ (seizures) [8]; the synchronization of neuronal activity in networks has therefore received significant attention, and network topology is considered to be a dominant factor affecting spreading dynamics [9–14], independent of the specific model by which activation is transmitted across the network. A wide range of transmission models has been employed in the literature. Representative examples include simple cellular automata-type spreading rules [11,15] and pulse-coupled oscillators (see [9] and references therein), of which integrate-and-fire models [13,16,17] are a special case. We note that in previous work, e.g. [9–11] (and in the present work), network nodes are to be interpreted as ‘neural units’ comprising many neurons, and representing a single cortical column, say. Detailed consideration of synaptic signalling models is therefore not appropriate.

A further essential aspect of neural systems is that they are spatially embedded—i.e. their nodes and edges are constrained to lie on a fixed geometric structure, and it is expected that these additional factors will further influence network dynamics. Indeed, though spatial networks have received considerable attention of late [18], the influence of spatial embedding upon network dynamics is not well understood.

A small number of recent studies have begun to employ neuroimaging data (e.g. via diffusion tensor imaging) to infer large-scale network connectivity [19,20]; however, the majority of studies in macroscale epilepsy and seizure modelling typically employs uniform lattices in one dimension or two dimensions with network connectivity restricted to nearest neighbours, or with additional long-range connections obeying arbitrary rules [17,21,22]. Such an approach fails to exploit a wealth of neuroimaging data [23], which reveals the intricate connectivity structure and surface geometry of the brain. Therefore, current approaches still observe simplified two-dimensional surfaces [24–26] inspired by the original studies of excitable media [27]—an organization that is different from the rounded shape of model organisms such as rodents or convoluted brain surfaces, such as for humans. Note that this limitation not only relates to epilepsy modelling but also to studies of spreading depression [28,29].

In this work, we seek to address the deficiency in neural network studies discussed earlier. We employ detailed spatial information obtained from freely available neuroimaging data of a rat brain to define a physiologically relevant neural network, and via comparison with commonly employed network architectures, determine the influence of the connectivity and complex surface geometry of the brain on seizure dynamics. We restrict attention to lattice-like network structures (i.e. networks that display significant clustering and long average path lengths), and a basic spreading model in favour of the more complex signalling models outlined above, as the aim here is to highlight the contribution of spatial network properties to signal transmission. We also limit our current study to spreading on the brain surface; the role of long-distance white matter fibre tracts between brain regions is not part of this work.

Our investigations highlight clearly that employing cortical surface geometry to inform network structure influences dramatically the propagation of activation through the network. By doing so, we indicate that activation dynamics of relevance to epileptic seizure initiation and progression (in particular, those highlighting the importance of the site of initiation within the network) display distinct differences in idealized networks, which are commonly employed in the literature. Most strikingly, we show that our cortical network delays significantly the total activation of the network when compared with idealized networks with the same coarse statistics: in the parameter regimes that we study, the time to activation is increased by a delay factor lying in the range 1.45–1.88. In this way, we highlight clearly the importance of realistic network structure on activation dynamics and indicate that such considerations should be included in theoretical models which aim to provide a more complete description, for example, of epileptic seizure activity.

The remainder of this paper is organized as follows. In §2, we use neuroimaging data to construct a lattice-like network architecture embedded on the cortical surface of a rat brain, together with a simple cellular automaton-like rule governing network activity. In §3, we present a comparative analysis of spreading dynamics over the cortical network, a uniform square lattice and an ensemble of two-dimensional geometric random graphs. A summary of our main results, together with a discussion of future avenues of investigation, is provided in §4.

## 2. Spatial complex networks

In this paper, we investigate the influence of network structure on the dynamics of the processes occurring on the networks, with specific application to neural signalling. We achieve this by comparing the spreading dynamics of simulated neural activation within a rat cortex with those obtained on commonly employed spatial networks. All of the networks examined herein are unweighted, undirected and without loops.

A plethora of tools and techniques for characterizing complex networks exists [30]; however, in the context of spatial networks not all of these measures remain relevant. An important feature of neural networks is the propensity for nearby nodes to connect, and thus, such graphs tend to exhibit high levels of local clustering. In addition, vertex reachability (defined as the ability to travel between nodes *i* and *j* following connections within the network) impacts significantly on spreading dynamics. In view of this, systematic comparison of the network models considered here will be effected via the clustering coefficient, which is given mathematically as
2.1where *C _{i}* denotes the probability that any two neighbours of node

*i*are connected and

*N*the number of nodes. Additionally, we employ the characteristic path length,

*L*, which is defined as the number of edges in the shortest path between two vertices, averaged over all pairs of vertices.

### 2.1. Rat cortical network

Spatial coordinates defining the cortical surface of the rat were obtained from the Caret software package [31] and processed using the Caret Matlab toolbox. Figure 1*a* shows the cortical surface of the left hemisphere of the rat brain, with typical neural activity spreading obtained via numerical simulation (see §3) superimposed.

Restricting to the left hemisphere, we construct a cortical rat network, with nodes positioned on the *N*_{rat} = 9623 available data points (figure 1*a*), by connecting nearest neighbours according to the following process: (i) a minimally connected nearest-neighbour network is defined via the triangulation provided by the Caret software package; (ii) vertex pairs are connected if they lie within a Euclidean distance *r* of each other; and (iii) connections are removed if the shortest path, calculated using the nearest-neighbour network defined in (i), between connected nodes exceeds a predefined number of steps, which is defined experimentally. While Euclidean distance is a reasonable estimate of physiological connection length in general [23], the third step is necessary to remove spurious edges which arise for large *r* that are near in the ambient space, yet distant as measured on the cortical surface: shortest path length, defined via a simple mesh triangulation provides a convenient method with which to measure this disparity.

### 2.2. Standard network models

To discern the importance of network structure on activation spreading dynamics, we compare the spreading dynamics observed on the cortical network defined in §2.1 with that observed on a uniform square lattice graph and on an ensemble of two-dimensional geometric random graphs [32], network structures typical of those employed in the literature [17,21,22].

To construct a random geometric graph, we place uniformly and independently *N* nodes at random on the unit square, and form connections between pairs of nodes according to Euclidean distance. Similarly, a square lattice graph is obtained by forming distance-dependent connections within an *N*-point uniform square lattice on the unit square. Figure 1*b*,*c* illustrates these networks and a typical simulation of neural activity.

For comparability with the cortical network outlined in §2.1, we choose *N* = 9604 (the closest square number to *N*_{rat}) and control the number of edges in the graph via , the Euclidean distance scaled on the surface area of the rat cortex. The surface area is calculated via the Caret Matlab toolbox to be *S* ≈ 221 (arbitrary units), which provides a scaling: Additionally, we impose periodic boundary conditions.

Figure 2 shows how the network measures *C* and *L*, with which we quantify differences between the three networks, vary with the connectivity distance, *r*. Ensemble measures for the random graph are calculated from 10 000 realizations. These results indicate that the networks exhibit distinctly different features when restricted to short-range connections (small *r*), while highly connected networks are comparable. It is noteworthy that the value of the mean clustering coefficient for the random graphs tends to the theoretical value
2.2

obtained in [33], where *Γ*(*x*) denotes the gamma function.

We remark that for *r* < 0.2 we find , and therefore no connections exist in the square lattice graph (figure 2*a*). At *r* = 0.2, the connectivity of the lattice graph corresponds to a four-point nearest-neighbour stencil (so that the connectivity matrix is analogous to a discrete Laplacian operator), and since nearest neighbours of a vertex are not neighbours of each other, it follows that *C* = 0. In addition, note that due to interplay between the regular discrete structure of the lattice and our rule based upon Euclidean distance with which to add new edges, variation of the clustering coefficient for the lattice graph is non-monotonic with respect to *r*. The mean degree of the other networks is *d* ≈ 5 (data not shown). Finally, we note that for *r* < 0.35, the random graphs are not fully connected, hence the absence of data in figure 2*b*.

## 3. Spreading dynamics

Propagation of activation within each neural network defined in §§2.1 and 2.2 was governed by a basic spreading model [11,15], summarized as follows.

Nodes *i* are restricted to exist in one of two states, *x*_{i}: active (*x _{i}* = 1) or inactive (

*x*= 0). Starting from an initial activation state, simulation operated in discrete timesteps; from one timestep to the next, an inactive node became activated (or an active node remained in the active state) if it was connected to at least

_{i}*m*active nodes. Initial conditions comprised a small region of activation (1% of the total nodes in the network) surrounding a node selected at random; the propagation of this activation through the network under our simple assumptions on spreading dynamics provides a convenient and compelling method with which to highlight the differences imparted by the networks under consideration. We choose the mean fraction of activated nodes as our key metric with which to investigate the different networks. Ensemble measures of network dynamics on the random graphs and rat network were constructed from 10 000 simulations; behaviour in the uniform lattice is identical for all initial activation positions.

The value of *m* places a lower bound on the connectivity of the network for which network activation can occur, and influences the speed of spreading of activation in the network (and, together with the value of *r*, the shape of the advancing activation front). As we consider highly simplified dynamics in this study, omitting, for example, random inactivation or complex intra-node dynamics (the better to emphasize the importance of network structure on activation dynamics), the balance between *m* and *r* determines completely the speed of activation of the network (indeed, for appropriate *m* and *r*, whole network excitation is inevitable) and, furthermore, affects all networks in the same manner. In all of the simulations that follow, we chose *m* = 2 without loss of generality. For *r* < 0.35, the network structures are not comparable (as highlighted by figure 2*b*); for *r* > 0.5, the networks are highly connected and activation spreads rapidly over the cortical surface. Our interest here is in the short-range connections between cortical columns, rather than long-range white matter fibre tracts. We therefore restrict attention to connection distances *r* ∈ [0.35, 0.5], leading to mean degrees in the cortical network 〈 *d*_{rat} 〉 ∈ [17.59, 37.36]. While actual estimates for mean degree in the rat cortex are not available, such connectivity is typical of that employed in the literature (see [9–11] and references therein), particularly when one accounts for the large spread of observed degrees (e.g. for *r =* 0.35, the degree ranges from 8 to 48 while for *r =* 0.5, it lies between 20 and 87), and so serves to illustrate our methodology. For brevity, in the figures that follow, we illustrate the differences in spreading dynamics obtained in each network for the choices *r* = 0.35 and *r =* 0.5.

Figure 1 shows a typical pattern of activation at an illustrative timepoint in each of the three networks, which serves to highlight clearly how differences in the underlying network connectivity are made manifest in the spreading of activation through the network. We remark that patterns of excitation are ‘well defined’ and uniform in the rat and lattice graphs, with all nodes within a region activated, while in the random graphs, lack of connectivity can lead to bottlenecks (here, we choose *r* = 0.25 to highlight these structural differences).

Figure 3 shows the spread of activation through each network, as measured by the fraction of activated nodes (averaged over all instances). Figure 3*a* indicates that for lower network connectivity (*r* = 0.35; the mean degree in each network was comparable: 〈*d*_{rat}〉 = 17.59, 〈*d*_{rand}〉 = 16.09 and *d*_{latt} = 20), the activation speed in each of the three networks differs significantly: the lattice graph requires dramatically fewer timesteps to achieve entire network activation; the cortical network is the slowest of the three. In fact, we observe a 1.88-fold and 1.45-fold increase in activation time in the rat network when compared with the lattice and random graphs, respectively. For more highly connected networks (*r* = 0.5), the cortical network remains the slowest to activate; however, the lattice and random graphs now display similar activation rates, with the uniform lattice being marginally faster. The delay imparted by the cortical network is now 1.65-fold (lattice) and 1.46-fold (random graphs).

Figures 4 and 5 highlight the influence of initial activation position on the spreading dynamics, which is quantified by the time to full-network activation, . For the values of connectivity distance, *r*, analysed here, full-network spreading was observed for all networks independent of topology. Therefore, all simulations contributed to the calculation of .

In figure 4*a–c*, histograms are presented, depicting the spread of observed in the simulations for each network; figure 4*d*–*f* shows the corresponding activation spread for each simulation, together with the ensemble mean. The delta function obtained in figure 4*b* and the corresponding results shown in figure 4*e* reflect the fact that the network dynamics are identical for all realizations owing to the uniform structure throughout the lattice and the periodic boundary conditions. By contrast, figure 4*a*,*d* indicates a very wide spread of activation times, and with no clear distribution, in the cortical network. In the random graphs, displays small variation on each realization; a wider spread is exhibited for a less well-connected network (data not shown). As highlighted by figure 4*c*, this distribution is well characterized by a Gaussian with mean * μ* = 43.03 and variance

*σ*= 0.6901 (

*p*< 0.01). The observed distributions of spreading dynamics for the cortical and the random networks (figure 4

*a*and 4

*c*, respectively) were found to differ significantly (

*p*< 10

^{−}^{5}over the investigated parameter range) according to the Kolmogorov–Smirnov statistic.

In figure 5, heat maps are presented to highlight those initial activation sites in the cortical network which provide significantly higher full-network activation speeds. These results indicate that the differences in activation time shown in figure 4 are due to the geometric structure of the rat brain and the presence of curved and smooth regions in the rat cortical surface.

The results shown in figure 5 are for the highly connected network (*r* = 0.5); however, qualitatively similar results are obtained for the range of parameter values analysed here: although values of vary, in all cases, initial activation of the folded regions leads to (approx.) a 1.7-fold reduction in total activation time (data not shown).

We have presented simulation results which highlight the differing rates of network activation, from an initial state comprising a small region of activated nodes. Our model can therefore be thought of as broadly applicable to the initial stages of epileptic partial seizures, whereby spreading can initiate from a certain region, leading to increased activity patterns in larger parts of the brain (mechanisms that lead to the rise of initial activity, such as high-frequency oscillations through gap junctions [34] or to seizure termination are not considered). Additionally, the results that we have presented indicating the dependence of activation progression on initiation site are pertinent to observations in epilepsy, as not all brain regions have the same probability of being starting points for seizures. We remark, however, that the concept of a specific focal point of seizure origin from which seizure activity spreads (owing to local abnormal connectivity) has been criticized [35], and alternative hypotheses presented [36,37]. We do not discuss this further as our focus here is on the global network structure imparted by the rat cortex embedding and its effect on activation propagation.

We remark that although our model has relevance to epileptic seizure spreading processes as described earlier, its form is intentionally simplistic: we have omitted a plethora of physiologically important network and signalling features (such as long-range connectivity or complex node activation dynamics) as our study addresses a fundamental aspect of the theoretical modelling of neural networks. Our central result is that cortical surface geometry affects profoundly signal propagation through the network, when compared with idealized networks with the same coarse statistics: the results contained in this section highlight clearly that significant differences in spreading dynamics are obtained in the network based on the cortical structure of the rat brain, compared with more standard network models. Since networks analogous to these standard models (defined in §2.2) are employed frequently in the theoretical literature in order to study the features of epileptic seizure initiation and dynamics, the simulation results highlighted by figures 3–5 indicate clearly that certain dynamics of relevance to seizure initiation and progression (in particular, those highlighting the importance of the site of initiation within the network) are not well captured by such approaches.

## 4. Discussion

In this paper, we have investigated numerically the influence of the structure of a spatially embedded complex network on the dynamics of the processes which occur upon it, with application to the development of improved theoretical models of the progression of epileptic seizures and spreading depression over the cortical surface.

The key feature of this work is that we employ neuroimaging data from the left hemisphere of a rat brain to define a neural network, whose spatial embedding represents accurately the structure of the cortical surface. Typically, theoretical studies of neural network dynamics employ uniform lattices or random graphs in one dimension or two dimensions: their focus being on in-depth analysis of various theoretical connectivity rules (e.g. small-world or scale-free networks) or signalling processes on the network dynamics. Here, we define connectivity within our cortical network by a simple criterion based on Euclidean distance (modified to account for the folded cortical structure), and employ a basic spreading rule to govern the propagation of activation. In this way, we highlight clearly the influence of physiologically relevant network structure on seizure dynamics, in isolation.

We compared numerical simulation of the propagation of neural activity within three different network architectures: the cortical neural network, a uniform square lattice graph and an ensemble of two-dimensional geometric random graphs, and studied these over a range of network connectivities. Note that we limit our current study on comparisons of spreading on the brain surface to short-range connections between nearby cortical columns; the role of long-distance fibre tracts between brain regions is not part of this work. In the parameter range chosen for the dynamic simulations which we investigate in detail, these networks were shown to be comparable in terms of their clustering coefficient and characteristic path length; in the case of networks with low connectivity (i.e. for values of *r* approaching the uniform lattice spacing), the networks display distinctly different characteristics. Despite their comparability, however, our simulations indicate dramatic differences in the propagation of activation through the various networks: in particular, we observe that the time taken to full activation in the cortical network is increased by a factor lying in the range 1.45–1.88. Especially striking is the importance of the site of activity initiation: a dramatic spread in the number of timesteps required to achieve full-network activation is observed in the rat network, even in the case of a highly connected network (*r* = 0.5, 〈*d*_{rat}〉 = 37.36). Small variation is observed in the random graph; the dynamics on the lattice graph is independent of initiation site. For less well-connected random graphs (*r* = 0.35, 〈*d*_{rand}〉 = 16.09), the variation is increased; however, the differences between all three networks remain significant (data not included). Moreover, our simulations highlight that, across the parameter range investigated here, different initial activation sites lead to an approximately 1.7-fold reduction in full-network activation time owing to variation in connectivity across these regions.

We employ a highly simplified model with which to investigate the spreading of activation over our network. Various features of relevance to physiological neural networks, such as complex node dynamics, long-range connectivity or stochasticity, are omitted from our formulation; indeed, our model draws no distinction between normal activation spreading and that seen in epilepsy. However, such an approach allows us to provide a more powerful exposition of the importance of the network structure imparted by the cortical geometry, in isolation. We have shown that certain features of the activation dynamics display distinct differences, when compared with idealized models of the type commonly employed in the literature. Most strikingly, we highlight a variability in activation time for the cortical network, depending on initial activation site. This is pertinent to observations in epilepsy, as not all brain regions have the same probability of being starting points for seizures. Variability in cortical networks is well studied in terms of structural and functional connectivity of brain regions, for example, concerning the degree of nodes [38–40] and its consequence on network robustness and performance [41,42]. However, our result is remarkable in that variability is observed for curved brain surfaces even in the absence of white matter fibre tract connectivity. Our results lead us to conclude that studies which do not take into account the spatial embedding of the cortex risk simplifying neural activation dynamics in a potentially significant way and are, therefore, unlikely to be able to represent accurately activation dynamics of relevance to epileptic seizures. We note, however, that in a physiological setting, similar disparity in network activation may be induced by including a range of other factors since the influences on the network dynamics within physiological neural networks are myriad; indeed, there is general agreement that no single factor can explain the varied phenomena associated with epileptic seizure dynamics.

We remark that epileptic seizure events are extremely rare in wild-type rodents, and characteristics of epileptic activity in animal models can differ from those observed in humans [34]. However, our initial study employing a network based upon the cortex of a healthy rat has highlighted the potential importance of geometric structure in activation dynamics; similar investigations in a network whose spatial embedding represents the convoluted human cortical surface therefore forms important ongoing work. In addition, future work will investigate how our predictions are altered under (i) a signal transmission model which represents more accurately the behaviour a neural unit and (ii) the introduction of long-range connections representing white matter structural connectivity.

## Acknowledgements

M.K. was supported by the WCU program through the KOSEF funded by the MEST (R31-10089), EPSRC (EP/G03950X/1) and the CARMEN e-science project (http://www.carmen.org.uk) funded by EPSRC (EP/E002331/1).

- Received January 8, 2013.
- Accepted January 24, 2013.

- © 2013 The Author(s) Published by the Royal Society. All rights reserved.