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).
 |
 |
,
, and 
inputs per msec. Lighter shades
denote higher discharge rates (maximum rate 120 Hz). Note the nearly
hexagonal pattern of clusters or ``hot spots" of activity. On the
right, we illustrate the motion of a typical cluster. Each vertex in
the graph represents a tracked cluster's position averaged over 50
msec. The diffusion constant is 
lattice units/msec.
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.
