The model consists of a cyclic lattice of units characterized by a voltage variable . Each unit models a nerve cell as a simple RC circuit ( msec) with the addition of a reset mechanism. Once a unit reaches a threshold voltage , it emits a pulse that is transmitted in one iteration (1 msec) to connected neighboring units, and the potential is reset by subtracting the voltage threshold. The unit dynamics are

where the current

represents the coupling of the unit to nearby unit through excitatory connections (matrix ) and inhibitory connections (matrix ). Additional input is provided by an independent external current . Such a model, composed of the ``integrate-and-fire'' [24] units with the dynamics described in eq. 1, is functionally similar to a recently studied model of self-organized criticality for earthquakes [38].

**Figure 1:** *
Center surround connectivity pattern. The cell
(not shown) at the center of the rectangular array is connected in a
probabilistic manner to units within a given distance determined by a
Gaussian distribution with lattice constants. These
short-range connections are excitatory (squares). The center cell
also inhibits a fixed fraction of cells on an annulus 8 and 9 lattice
constants away (triangles). During a particular simulation, the
connectivity pattern is fixed, although the exact synaptic weight
varies stochastically. *

Based on physiological evidence, the pattern of local excitatory and
inhibitory connections is * center-surround*, for the validity of
this assumption can be found. Hess and coworkers [19]
found that iontophoretic application of the excitatory agonist
glutamate to visual cortex of anesthetized cats induces excitation of
neurons within of the application site and distant
inhibition at distances between 100 and . Similar studies
in rat somatosensory cortical slice preparations
[43] confirm the general pattern of an inhibitory
surround enclosing the excited region. Each unit is excitatorily
connected to **N = 50** units chosen from a Gaussian probability
distribution of lattice constants, centered at the
unit's position (see Fig. 1). Likewise, **N**
inhibitory connections per unit are chosen from a uniform probability
distribution on a ring eight to nine lattice constants away. The
weight of the excitation and inhibition, in units of voltage
threshold, is and . External
input is modeled independently for each cell as a Poisson process of
excitatory pulses of magnitude , arriving at a mean rate
.

Scaling the relative degree of inhibition while keeping the sum of excitation and inhibition constant leads to a transition from a spatially homogeneous state to a clustered activity state. This transition can be understood using a mean-field description of the dynamics, where we write the pulse rate as a function of the averaged currents : [2], where is the minimum dead time between pulses. In this approximation, the dynamics associated with eq. 1 simplify to

The homogeneous solution is stable only if the connectivity pattern
satisfies for all **k**, where
is the Fourier transform of . As
one increases the relative strength of inhibition, clusters of high
firing activity develop. The clusters form a hexagonal grid across
the network, as illustrated in Fig. 2a. Clusters
do not remain stationary but pursue a stochastic motion over the layer
(Fig. 2b).

The level of external input also influences the spatial pattern; for even higher values of input the clusters merge into stripes. The transition from a homogeneous state to hexagonal clusters to stripes is generic to many nonequilibrium systems in fluid mechanics, nonlinear optics, reaction-diffusion systems, and biology [8,9,12,36]. In the mean-field description, the dynamics (3) are subject to a Lyapunov potential function [21] so that the system must relax to a stable state---persistent dynamics are ruled out.

However, the mean-field equations
[1,52,53] do * not*
capture the essential dynamics of pulse-coupled networks responsible
for the spatial pattern's metastability. Partial synchronization
inside local populations induce large fluctuations which destabilize
the mean-field pattern. Indeed, simulations performed with constant
(instead of stochastic) input and with isotropic connectivity reveal
that the system is intrinsically chaotic. These fluctuations allow
clusters of high activity to diffuse throughout the system, as can be
seen from the cluster's trajectory in Fig. 2b.
Repulsive interactions with surrounding clusters generally constrain
the motion of a cluster to remain within a certain radius. This
vibratory motion is occasionally punctuated by longer-range diffusion.
If one tracks the activity of individual cells during a period of
time, one observes large fluctuations leading to high variability in
interspike intervals and in the total number of spikes, which are due
to the fact that even under stationary stimulation the system
generates its own rate fluctuations.

**Figure 3:** *
Coefficient of variation, , of a representative range
of cells shown against the average time between spikes, that is, the
inverse of the mean firing rate. The solid dots in the upper cloud are
from our model. The crosses in the middle cloud represent the
behavior when all network effects are eliminated (i.e. ),
and units receive both an excitatory and an inhibitory stream of
external Poisson input, with . Without
inhibition, the associated values are reduced by a factor of
about 2. The diamonds in the lowest cloud are from cells in a random
network with sparse, non-local connections that have no organized
topography. The same number of inhibitory and excitatory connections,
50, were used as in the local center-surround connection scheme. *

The coefficient of variation () computed for all the cells in the network, over a 400 sec long simulation are displayed in Fig. 3 Note that almost all values of are on the order of one or larger. The pattern of observed 's reproduces qualitatively the values measured for cells in cortical areas V1 and MT in the awake monkey responding to bars and to clouds of moving dots [44].

The coefficient of variation does not contain enough information,
however, to fully characterize the statistical properties of the spike
train. In order to obtain a more complete characterization of those
properties, we computed the inter-spike-interval (ISI) distribution
and the cross-correlation function of spike trains
(Fig. 4). We find that the ISI distribution decays as a
power law in **t** with exponent over one decade (from 25
to 300 msec). The correlation decays slowly via a power law
with exponent , as shown in Fig. 5. At the upper
temporal cutoff of 300 msec, the correlation function reaches the
baseline of chance coincidence. The power spectrum at low frequencies
has a form consistent with the correlation and decays as
(see, for more details [56]).

**Figure 4:** *
The average ISI distribution for 19 cells. Superposition of
exponential distributions from the cluster frame of reference (solid
line) results in a power law in the fixed unit
frame. *

**Figure 5:** *
Measured and predicted (solid line) correlation on a
log-log scale. The correlation function decays as due to
the slow diffusion of activity clusters.*

Such power-law distributions of the ISI and correlation functions have been reported in neurons from the somatosensory cortex [60], in the auditory pathway [49,50], in the mesencephalic reticular formation [18] and in visual cortex [48]. One should note that such power-law functions are not explained by the variability of spike-trains described by Poisson processes (which have exponentially distributed ISI) and which are predicted by models based on a precise balance of the excitatory and the inhibitory inputs to cortical cells [42,54]. Another characteristic property of the data which is explained by our cluster-diffusion model is the typical ``castle-on-hill'' form of cross-correlations: [33] found cross-correlations which show a narrow peak 5-50 msec, riding on a much wider, 100-1000 msec, ``hill''. In our model, the narrow peaks result from the synchronized activity of cells inside a cluster while the wide hills from the much slower diffusion of the clusters.

A much more direct method to test the prediction of our model is by means of imaging techniques that tap directly into the spatio-temporal patterns of neural activity. Using the novel technique of optical imaging with fast voltage sensitive dyes [39] report moving elliptical foci of neural activation in response to electrical stimulation. Moreover, in a recent study Arieli et al. (1995) [3] have shown that moving patterns of spontaneous activity on the cortical surface contribute to the variability of neural response.