Wave Speed in Excitable Random Networks with Spatially Constrained Connections

Very fast oscillations (VFO) in neocortex are widely observed before epileptic seizures, and there is growing evidence that they are caused by networks of pyramidal neurons connected by gap junctions between their axons. We are motivated by the spatio-temporal waves of activity recorded using electrocorticography (ECoG), and study the speed of activity propagation through a network of neurons axonally coupled by gap junctions. We simulate wave propagation by excitable cellular automata (CA) on random (Erdös-Rényi) networks of special type, with spatially constrained connections. From the cellular automaton model, we derive a mean field theory to predict wave propagation. The governing equation resolved by the Fisher-Kolmogorov PDE fails to describe wave speed. A new (hyperbolic) PDE is suggested, which provides adequate wave speed that saturates with network degree , in agreement with intuitive expectations and CA simulations. We further show that the maximum length of connection is a much better predictor of the wave speed than the mean length. When tested in networks with various degree distributions, wave speeds are found to strongly depend on the ratio of network moments rather than on mean degree , which is explained by general network theory. The wave speeds are strikingly similar in a diverse set of networks, including regular, Poisson, exponential and power law distributions, supporting our theory for various network topologies. Our results suggest practical predictions for networks of electrically coupled neurons, and our mean field method can be readily applied for a wide class of similar problems, such as spread of epidemics through spatial networks.


Introduction
Different types of networks are found across many scales, from metabolic networks in a single cell, to neural networks in brain, up to social and technological global networks. The theory of networks receives increasing attention since the pioneering works that formulated random graphs [1], and the recently discovered ubiquity of small-world networks [2] and scale-free networks [3]. Reviews on general theory of networks can be found in [4][5][6]. A comprehensive up-to-date review of spatial networks is given in [7].
Since its first formulation [1], the Erdös-Rényi (ER) graph became a cornerstone of network theory. An ER graph G NM consists of N nodes and M links (edges), and each link connects two nodes which are selected randomly. In a sufficiently large network, the number of links emanating from a node (degree) is a random variable with Poisson distribution P(k)~SkT k e {SkT =k!, where SkT is the network mean degree, SkT~2M=N. Despite its advantages, the ER graph is not suitable for studying spatial phenomena because it is spatially homogeneous. However, in most real-world networks the connections are spatial and variable in length. Also, the maximum length of connection is usually limited by the available resources or other natural restrictions. To address this problem, spatial generalizations of the ER graph were suggested. For example, two nodes can be connected only if the distance between them is below threshold r c [8]. This model was used to simulate spatio-temporal activity in networks of electrically coupled neurons [9]. Another example is the Waxman model, in which the probability that two nodes are connected is a decreasing function of distance between the nodes [10]. The latter model was used to simulate the Internet [11].
In many networks the nodes are excitable, meaning that active state can arise and propagate from one node to another if they are connected. In this way, action potentials propagate through neural networks, computer viruses spread in the Internet, and diseases are transmitted through transport networks. If the nodes are excitable, dynamical states propagate through a network both temporally and spatially, leading to waves and more complex patterns.
A case study in our work is the emergence of spatiotemporal patterns with very fast oscillations (VFO, w80 Hz) measured by electrocorticography [9], recorded in neocortex of patients prior to epileptic seizures ( Figure 1A). There is growing experimental and theoretical evidence that VFO are caused by electrically coupled pyramidal neurons which are connected by gap junctions, thus providing direct excitation from one to another, which does not require synaptic transmission [9,12,13].
Nodes in such a network are dynamical excitable units. Although their intrinsic behavior can be complex and require detailed multi-compartment model of each neuron [14][15][16], understanding of the neural network spatiotemporal oscillations requires highly reduced models, such as cellular automata which capture three main states (resting, firing, refractory) with a minimal set of parameters [13]. The system can be described by a network version of Greenberg-Hastings cellular automaton (GHCA), a discrete model of an excitable medium [17]. Variations of GHCA have been used in many studies, including collective oscillations of pyramidal cells in the hippocampus [8,13,18], sensory networks [19,20], and the evolution of HIV infection [21].
In the network version of GHCA that we use, cellular interactions can be distant rather than next-neighbor [8]. The cells are connected into Erdös-Rényi random graph with spatially constrained connections (hereafter SCC), so the distance between connected nodes is not greater than connectivity radius r c . Under random spontaneous activation of some of the resting cells, large networks demonstrate oscillatory dynamics with complex spatiotemporal activity driven by many interacting waves. Traub and coauthors [9] have recently demonstrated these patterns to be strikingly similar to those observed in ECoG recordings ( Figure 1B). In the model, the complex patterns of activity arise when multiple waves are born from spontaneously activated cells and they grow and coalesce in time and space.
A single active node may generate an expanding circular wave of excitation, if the mean network connectivity is high enough. The wave maintains its shape and travels with constant speed, which is an important characteristic of system's excitability and depends on network topology. Knowledge of wave speed in excitable networks allows prediction of how fast the active state (neuronal activity, viral infection) propagates and invades the rest of the network. Although simulations can be done in each particular case, it is important to have a basic mean field theory Figure 1. Neural network activity in experiments and in the cellular automaton model. A. A snapshot of electrocorticographic (ECoG) data of brain activity, measured by 866 subdural array of electrodes. Data is interpolated between nodes, white areas correspond to high activity. B. A snapshot of activity from a cellular automaton model in an 4006400 network. The network is subject to noisy input from spontaneously activating cells (rate 1:25|10 {5 ). Active cells are white, refractory and excitable are black (simplified color code). C. Snapshot of activity in a 10610 sub-network with detailed color code: red for active, blue for refractory, black for excitable nodes. Lines show links between nodes. D. Rules of the CA model: excitable node (black) may become active (red), if activated by a neighbor. After being activated, the node becomes refractory (blue) for a period of time t r , after which it becomes excitable again. Data in A is a courtesy of Miles Whittington, recorded in Patient B of [29]. doi:10.1371/journal.pone.0020536.g001 (MFT) to understand and predict the network dynamics without simulations.
In this paper we derive an MFT to predict wave speed in a random (Erdös-Rényi type) network with spatially constrained connections (SCC), for given connectivity radius r c and mean degree SkT. The results are generalized to Erdös-Rényi type networks with various radii distributions, and further to non-Erdös-Rényi networks with various degree distributions, suggesting universality of our mean field theory.

The system
In the following section we describe the system and its dynamic properties, along with simulations that demonstrate its typical behavior.
2D network. We study a 2D network of excitable cells (nodes), which are are set on a uniform N x |N y grid, with unit space between adjacent nodes ( Figure 1C). Connections (links) are bidirectional, i.e. activity can be transmitted in both directions. Links can be long but limited by r c : within a circle of radius r c (network with round 'footprint') or within a square with side r c (network with square 'footprint') [8]. Generally, the connectivity radius r c is much larger than 1, but much smaller than array dimensions N x (N y ). The number of links per node (degree) follows the Poisson distribution P(k)~SkT k e {SkT =k!, whereas the distribution of link lengths is uniform in (1, . . . ,r c ), unless stated otherwise.
Quasi-1D network. To treat the system analytically, we reduce it to a quasi one-dimensional network: the length of links is limited only along X but may be unlimited in the Y coordinate. In other words, it is a 2D network which has an interval footprint (Dx{x 0 Dƒr c ). A quasi-1D network provides the same wave speed along X as the 2D network with square footprint (with same r c ), which is also consistent with 2D network with round footprint (simulations not shown).
Cellular automaton model. A node in excitable state (E) becomes firing (F ) if one or more of its connected neighbors in the network are firing. After one time step Dt the firing node becomes refractory (R 1 , . . . ,R tr ) for a relatively long period t r , after which it becomes excitable again. Thus each node rests in (E) or undergoes a sequence of states (E,F ,R 1 , . . . ,R tr ,E) if activated by a neighbor ( Figure 1D). Initially all nodes are in excitable state, except a small number of firing nodes that initiate the wave. Node states are updated simultaneously at each time step.
Initiation of wave in a small 2D network is shown in Figure 2 (first four time steps). Although directions of links are random, activity propagates outwards from the initial point because it is followed by refractory state, prohibiting backward propagation. Propagation of waves in large 2D networks is shown in Figure 3A-C. Waves in networks with round and square footprint do not differ qualitatively (not shown), because neighbors of each node are chosen as random lattice points from round or square neighborhood (respectively), which differ only in relatively small corner areas.
In the rest of the paper, we work with quasi-1D networks, because they allow transparent mean field analysis, and yet behave almost identically to 2D networks, with the only difference that wave fronts are linear rather than circular ( Figure 3B vs. A). Next, we start from simpler case where all links have maximum length (r~r c ), and then analyse the case where links are distributed randomly (0vrƒr c ).

Mean field theory
Links of maximum length. In order to treat the system analytically, we assume that the network is regular (all nodes have same degree k), and all links have same length (r~r c ). Later we will show that these crude assumptions provide a good approximation for more general cases.
We assume that each excitable cell can become firing by one (and only one) of its neighbors that is firing, which is true for wave front where firing cells are rare. Each firing cell is surrounded by excitable cells, and may produce at maximum (k{1) firing neighbors. One of k links is missing because it points to a cell from which firing has come. In other words, if a cell is excitable (probability E), it can be activated by one of its neighbors (probability SF T, spatial average of F). Otherwise, if a cell is already firing (F ) and there are excitable neighbors (prob. SET), it produces (k{1) firing cells. These relations can be written as where F and E are the probabilities of a cell to be in firing or excitable state, respectively, SF T and SET are spatially averaged F and E in the neighborhood of x. This equation holds only for wave front, where firing cells are scarce, excitable cells are abundant, and there are no refractory cells yet.
With the notion that most cells are excitable (F&0, E~1{F &1), the E-terms become unitary (linearized equation).
In the quasi-1D system, SF T~1 2 (F (x{r,t)zF (xzr,t)), and the equation simplifies to Taylor's expansion in x at the right hand side gives Taylor's expansion in t at the left hand side gives Parabolic (Fisher-Kolmogorov) equation. Taking into account only the first time derivative, we arrive at a linearized version of Fisher-Kolmogorov equation Dt is the growth rate, and the second-order extinction term is omitted for wave speed analysis. The well-known [23] formula for Fisher- gives infinite growth of speed with mean degree k ( Figure 4, upper line). Taking high-order terms into account in the right hand side of the PDE does not alter the principal behavior of wave speed (simulations not shown), demonstrating that parabolic PDEs are not suitable for wave speed prediction. Hyperbolic equation. Keeping both first and second time derivatives in the left part of Eqn. 4, we obtain a hyperbolic equation for firing node density Wave speed can be found by marginal stability analysis [24]. Substituting variable z~x{vt, we obtain the equation The roots of the characteristic equation must be real, which gives the minimum wave speed We put the time step Dt~1 as in the CA model, and k~SkT. This wave speed demonstrates qualitative agreement with the CA model on an ER SCC network (Figure 4, lower line). Most importantly, v(k) gradually saturates to the maximum possible speed v(k)?r c for high k&1, in agreement with CA simulations and intuitive expectations.
The wave speed v(k) is shown in more detail in Figure 5 (upper solid line). One can see that v(k) falls near to CA simulations on networks where all links have maximum length r c ( Figure 5, triangles), as expected. Surprisingly, v(k) also approximates well the CA simulations on networks where links have random length 0vrƒr c ( Figure 5, circles), which is our primary model. This phenomenon is explained in the next section.

MFT for links with random length
The case where links have random length 0vrƒr c is derived similarly to the maximum-length case shown above. Recall that the spatial scale between nodes is unity, r c &1 and r~1 . . . r c . Let all lengths be equally probable, P(r)~1=r c . All nodes have degree k, so the fraction of nodes with at least one link of length r is kP(r)~k r c . The evolution of firing nodes density in time is given by Integration over 0vrƒr c and omitting E*1 gives Taylor's expansion up to second derivative of the function under the integral gives the equation The latter equation is equivalent to Eqn. 3, with the only difference being a reduced radius (r~r c ffiffi ffi 3 p ). Therefore, the wave speed for the random-length links would be simply scaled by ffiffi ffi 3 p However, this naive scaling based on uniform radii distribution underestimates the wave speed of CA roughly by a factor of two ( Figure 5, lower solid line, compare to circles). This happens because cells in the wave front have actually non-uniform distribution of links from which they have received their activation ( Figure 6A). The distribution is strongly biased towards the longest links. A wave front generates a new wave front at the next time step by sending activity through the longest links out of available k. To support this notion, we generated sets of k i.i.d. discrete random variables r uniformly distributed in (1, . . . ,r c ) and computed their mean maxima Smax(r 1 , . . . ,r k )T as a plausible estimate of wave speed, that is a contribution of activity propagation from each single node to a global propagation of wave per unit time. As one can see in Figure 6B (broken line), the mean maxima of k random radii gives a good measure of CA wave speed (circles), especially for high k. These numerical calculations are supported by analytic formula for the expected mean of maxima (black solid line), Smax(r 1 , . . . ,r k )T~r c 2 {1=k , derived for a uniform continuous distribution r [ (0,r c ). This formula gives a very good prediction of CA wave speed at high k, demonstrating that the wave propagation is indeed mainly determined by the mean maximum of k radii, which converges in the limit (k??) to the maximum possible radius r c . In other words, it is not the average, but rather the maximum link length that determines the wave speed in a random network with random radii distribution.
Role of link length distribution. To study the effects of other possible radii distributions, we simulated networks with five distributions ( Figure 7A): the uniform radii distribution, the fixedvalue distribution (r:r c ), a bell-shaped and two exponential distributions (increasing and decreasing, respectively). The speeds of wave propagation in resulting random networks are qualitatively similar and always significantly higher than the average value of the corresponding distribution ( Figure 7B). Broken lines in Figure 7B are mean maxima of k radii generated from each distribution, which served as a good estimate of wave speed in four out of five distributions. Note that the distribution for which our maximum link hypothesis works the least is the one with an exponentially small probability of reaching the maximum link length r c (see the purple line in Figure 7A).
Although we can not provide an analytic estimate of the effective radius in each distribution, these simulations support the notion that wave speed is predominantly determined by the maximum rather than average radius, with maximum taken out of SkT radius realizations.

High-order analysis
The use of Taylor's approximation in the PDE derivation could potentially bring unwanted errors because Dt and r are in fact not small. We use a high-order analysis to estimate wave speed in exact terms, by looking for a traveling wave solution in the form F (x,t)~e {l(x{vt) without using Taylor's expansion. See Methods and Algorithms for details. Figure 5 shows that wave speed v Ã (upper dashed line) obtained by this method is close to v(k) obtained from a hyperbolic PDE (upper solid line). This confirms that considering derivatives of order above 2 (in time and space) will not affect the wave speed qualitatively. Therefore, the hyperbolic PDE is both necessary and sufficient for qualitative wave speed prediction.

Role of degree distribution
The variation of node degrees depends on degree distribution of a network and strongly affects the wave speed. To ultimately simplify the network and add the variance in stages, we simulated the CA model on spatial networks the following degree distributions: N regular network (each node has same degree k), N fk{1,k,kz1g distribution: a node's degree is chosen randomly from fk{1,k,kz1g,k §2, In all networks, the lengths of links between nodes were uniformly distributed in (0vrƒr c ). As one can see in Figure 8A, the wave speed profiles vary widely when plotted against network mean degree SkT. However, they merge into nearly the same shape when plotted against Sk 2 T=SkT ( Figure 8B, inset). The key role of Sk 2 T=SkT ratio is evident from general network theory. For a randomly chosen link, the degree of a node on its end follows the nearest-neighbor distribution P nn (k)~k SkT P(k), where P(k) is the original network degree distribution [6]. The mean degree of a connected node SkT nn is therefore different from a randomly picked node: SkT nn~P kP nn (k)~S k 2 T SkT , so our ratio is merely the nearest-neighbor mean degree. In other words, when activation travels from one node to another, the degree of a node in the end of a link follows the nearest-neighbor distribution P nn (k), which has a mean of SkT nn~S k 2 T SkT and gives the actual branching ratio of the activation in the new node. For an ER graph, SkT nn~S kTz1. Substitution of SkTS kT nn {1 into wave speed (7) gives us a more general formula for fitting wave speeds  where SkT nn~S k 2 T SkT as before. As shown in (Figure 8, inset), this formula gives a good estimate of wave speeds in a variety of network degree distributions.

Discussion
In summary, we have analysed behavior of excitable random networks with spatially constrained connections (SCC), and derived a mean field theory of the activity propagation. We conclude that the hyperbolic PDE is necessary and sufficient to capture the wave speed in random (uncorrelated) networks with spatially constrained connections. The wave speed is mainly determined by the longest possible connection between the firing node and a node it can activate, so the mean maximum of SkT random radii Smax(r 1 , . . . ,r k )T serves as a good predictor of wave speed in Erdös-Rényi SCC networks with various radii distributions.
We have derived formula (7) for wave speed v(SkT), which agrees with simulated behavior of CA on Erdös-Rényi SCC networks. Simulations of CA on networks with other (non-Poisson) degree distributions suggest a more general formula (11) for wave speed, which depends on SkT nn~S k 2 T SkT , the nearest-neighbor mean degree. So, our mean field theory extends to networks with various degree distributions, provided that Sk 2 T SkT is used as a more universal measure of network's average branching, rather than simple SkT, which is explained by general network theory. The original [8,9,13] cellular automaton model of very fast brain oscillations assumed that pyramidal neurons are connected via gap junctions into an Erdös-Rényi SCC network. Our results show that wave speed (hence network excitability) will scale with radius r c and degree SkT nn the same way even if the network  topology is different. This has experimental implications: accounting for the wave speed does not require reconstruction of full gap junction network topology. Rather, determining the mean SkT and variance (s 2~S k 2 T{SkT 2 ) of the gap junctions per cell, and the maximum inter-neuronal distance (r c ) for coupled cells should suffice. Also, our simulations show that wave speed is primarily driven by the longest of SkT connections emanating from an average node. Therefore knowledge of the longest axoaxonal connections is more important than recovering their full range.
Our present knowledge of wave speed will facilitate a detailed analysis of waves interaction and coalescence, which generate complex oscillatory dynamics in large networks shown in [9].
In general perspective, the fact that wave speed depends on Sk 2 T SkT shows that variation in node degree plays an important role in wave propagation. High variation of node degrees provides high wave speed, presumably due to presence of highly-connected nodes (hubs). The hyperbolic PDE (12) should not be confused with similarlooking telegrapher's equation, because the (k{1) term in Eqn 12 is positive and gives a self-sustaining wave in time, in contrast to the decaying solution of telegrapher's equation.
The special type of networks we study (with spatially constrained connections) should not be confused with small-world networks. The limited length of link is crucial for spatial phenomena, whereas in conventional small-world networks the ''short-cuts'' are unlimited in length. This makes wave speed infinite and destroys spatial coherence, but improves temporal coherence [25,26]. However, small-world networks might be constructed with spatially constrained shortcuts, which could be an interesting system for analysis.
The CA model we use is similar to epidemiological SIRS model (susceptible-infected-recovery-susceptible). Therefore, our theory may help predict spatial spread of epidemics in large networks. For example, it may be applied to predict spatial spread of SIR type malware through a sufficiently large network of WiFi routers [27], where the length of connections is limited by the router's range.
We studied networks with degree distributions that cover a broad spectrum of possible Sk 2 T versus SkT combinations ( Figure 8B, inset), suggesting that our MFT holds for uncorrelated networks of arbitrary degree distribution (with spatially constrained connections). The degree distributions studied here appear in real-world networks, such as the neuronal network of C. elegans, the power grid network, acquaintance networks and the WWW (see [4] and references therein). Our MFT apply to these and other networks, provided there is a metric to measure connection length, and connections are limited by some constant r c .

Network with spatially constrained connections
The network consists of excitable nodes, which are are set on a uniform N x |N y grid, with unit space between adjacent nodes. The network with defined degree distribution P(k) and connectivity radius r c is constructed by a procedure similar to that for spatially homogeneous networks [28]. Initially, each node is assigned to a random number of 'stubs' k, which is picked from distribution P(k). Next, the program picks nodes from a randomized (shuffled) list of all nodes. The list of nodes must be randomized to avoid artificial correlations imposed by node order. For each picked node with nonzero number of stubs, the program randomly picks one of its neighbor within distance r c . If the neighbor has nonzero number of stubs, too, both nodes are linked. Their numbers of stubs are decremented by -1, and their numbers of links are incremented by +1. The procedure is repeated until all stubs of the chosen node become links. To avoid infinitely long search in a situation when all neighbors are already linked and they have no more free stubs, the search of potential neighbors is stopped after 10 4 unsuccessful attempts. The procedure is repeated for each node.
Single networks of size up to 400|400 were simulated in Matlab. Large-scale simulations of multiple (64 to 128) networks of size 1000|1000 with r c~2 0 (used in Figure 8) were carried out in C/C++ MPI program on an IBM Blue Gene supercomputer.

High-order analysis
The use of Taylor's approximation in space and time could potentially bring unwanted errors because Dt and r are in fact not small. Here we estimate wave speed in exact terms by looking for a traveling wave solution in the form F (x,t)~e {l(x{vt) without using Taylor's expansion, thus taking all high-order terms into account. Both cases of links with maximum and random length are analyzed below.
Maximum length. and can be solved numerically for marginal stability analysis. In order to find the minimal wave speed, we need to find v Ã such that the function in the left part f 1 (l r )~e lrv Ã has a unique common point with the function in the right part f 2 (l r )1 2 (e {lr ze lr )z(k{1), for given k. This happens when f 1~f2 and f 1'~f2' (plots touch each other). Parameter v Ã is the minimal wave speed, normalized to the radius r. Figure 5 shows that v Ã (upper dashed line) obtained by this method is close to v(k) obtained from a hyperbolic PDE (upper solid line), so considering derivatives of order above 2 (in time and space) will not affect the wave speed qualitatively. Therefore, the hyperbolic PDE is both necessary and sufficient for qualitative wave speed prediction. Random length. This case is analyzed in a similar way. Eqn.
This equation is solved numerically in the same way as in the maximum-length case. As one can see in Figure 5, the obtained minimal wave speed v Ã (lower dashed line) is close to v(k) obtained in hyperbolic PDE (lower solid line). However, as shown in Results, the MFT for maximum link length provides better approximation of CA simulations, so this case is shown here only for completeness of analysis.

Numerical integration of PDE system
Analytic formulas for wave speeds are supported by numerical integration of corresponding PDE systems. In the hyperbolic case, the closed nonlinear system reads By the normalization F zRzE~1 the equation for E is not necessary. Parameters are Dt~1,r~1,D~r 2 2Dt ,t r~1 5. This system was replaced by a system of three parabolic PDEs (by introducing G~F t ), which was solved by a method of lines in Matlab, with parameters dx~0:1, dt~0:002, mesh size 200|200, and explicit Euler scheme. Initial conditions were F~R~0 in all points except the center line where F~1.
Simulations of traveling waves in the Fisher-Kolmogorov equation were carried out similarly, with the first equation changed to F t~D F xx z(k{1)F (1{F{R).
Both systems demonstrate the effect of wave coalescence -two waves cancel each other at the areas where they meet, while their outer borders merge into one big circular wave. The main difference between the two PDE systems is the speed of wave, which grows infinitely with k in the Fisher-Kolmogorov system, but saturates to unity for the hyperbolic PDE system, as discussed in the Results.