## Abstract

A neuronal network inspired by the anatomy of the cerebral cortex was simulated to study the self-organization of spiking neurons into neuronal groups. The network consisted of 100 000 reentrantly interconnected neurons exhibiting known types of cortical firing patterns, receptor kinetics, short-term plasticity and long-term spike-timing-dependent plasticity (STDP), as well as a distribution of axonal conduction delays. The dynamics of the network allowed us to study the fine temporal structure of emerging firing patterns with millisecond resolution. We found that the interplay between STDP and conduction delays gave rise to the spontaneous formation of neuronal groups — sets of strongly connected neurons capable of firing time-locked, although not necessarily synchronous, spikes. Despite the noise present in the model, such groups repeatedly generated patterns of activity with millisecond spike-timing precision. Exploration of the model allowed us to characterize various group properties, including spatial distribution, size, growth, rate of birth, lifespan, and persistence in the presence of synaptic turnover. Localized coherent input resulted in shifts of receptive and projective fields in the model similar to those observed *in vivo*.

## Introduction

The cerebral cortex includes a population of billions of neurons, each making thousands of synaptic contacts with its neighbors. Given the complexity of the connectivity inherent in cortical anatomy, efforts to describe the pattern of electrical activity in exact detail within even a highly localized population of cortical neurons would be extremely difficult. Our understanding of cortical dynamics is also rendered difficult by the temporal properties intrinsic to such a biological system. For example, different types of neocortical neurons can generate a variety of spiking patterns in response to the same input (Connors and Gutnick, 1990; Gibson *et al.*, 1999). In addition, spikes are sent along neurites at variable velocities that give rise to axonal conduction delays that can be as large as 40 ms (Waxman and Bennett, 1972; Swadlow, 1991, 1992, 1994; Salami *et al.*, 2003). Moreover, empirical data indicate that the efficacy of a synapse can change rapidly in response to the precise timing of firing of pre- and postsynaptic neurons, a phenomenon referred to as spike-timing-dependent plasticity (STDP) (see Markram *et al.*, 1997; Bi and Poo, 1998; Debanne *et al.*, 1998; Feldman, 2000; Sjostrom *et al.*, 2001; Froemke and Dan, 2002; Izhikevich and Desai, 2003). As a result, such a system can exhibit quite rich spatiotemporal patterns of activity with a possibility that neurons spontaneously self-organize into neuronal groups (Edelman, 1987, 1993).

In order to gain some insight into the dynamic patterns that might be engendered by a network constructed according to the general architectural principles of neocortex, Pearson *et al.* (1987) simulated a model of a cortical somatosensory area. They showed that patterned input resulted in the formation of neuronal groups within the simulated cortical sheet, in accord with the theory of neuronal group selection (TNGS; Edelman, 1987, 1993). The simulation by Pearson *et al.* (1987) did not, however, explore the detailed dynamics of groups of spiking neurons with millisecond spike-timing precision.

In the present study, we describe a more detailed cortical simulation composed of neurons exhibiting the known types of spiking behavior found in mammalian neocortex. The model simulates the activity of 100 000 neurons and the changes in synaptic strength of 8.5 million synaptic contacts with sub-millisecond time resolution. As in the neocortex, the simulated neurons interact via both local and long-distance connections, and consequently generate integrated patterns of spiking activity via reentrant signaling (Edelman, 1987, 1993). Our goal was to study the fine spatiotemporal structures of emerging firing patterns. We found that, if axonal conduction delays and STDP were incorporated into the model, neurons in the network spontaneously self-organized into neuronal groups even in the absence of correlated input. Each such group is made up of tens to hundreds of neurons that can fire time-locked spiking patterns with millisecond precision.

Neurons in the model did not fire unless the network received some level of noisy input. While low input levels resulted in rhythmic synchronized bursts of spatially uniform activity, higher input levels gave rise to asynchronous Poissonian patterns of spiking activity that resembled those recorded *in vivo*. Although initially no regularities were discerned in the record of spike timing in this state, within a short period of time, certain recurrent firing patterns appeared spontaneously among subsets of neurons that were located close to one another. Spiking activity within these subsets was not synchronous, but their spikes were time-locked to one another, and thus they constituted spike-timing neuronal groups. Precise firing sequences of such neuronal groups resembled those recorded *in vivo* (Chang *et al.*, 2000; Lindsey *et al.*, 1997; Mao *et al.*, 2001; Prut *et al.*, 1998; Villa *et al.*, 1999; Tetko and Villa, 2001a,b). The neuronal groups appeared and persisted for variable periods of time as their constituent neurons were linked by strong synaptic connections with matching conduction delays.

By taking account of the detailed properties of the network, we were able to characterize properties of such neuronal groups including their spatial distribution, size, growth, rate of birth and lifespan. In the light of observations on shifting receptive field boundaries within somatosensory maps (Merzenich *et al.*, 1983), we also investigated how these group properties were altered when the simulated network received localized coherent excitatory input to simulate the effect of coherent sensory cortical stimulation.

## Materials and Methods

### Anatomy

As compared to real cortices, the model is obviously greatly reduced in the number of its neurons and synapses as well as in its anatomical complexity (see Fig. 1). Nevertheless, we made efforts to preserve important ratios and relative distances found in the mammalian cortex (Braitenberg and Schuz, 1991). The ratio of excitatory to inhibitory neurons was 4/1. The span of local non-myelinated axonal collaterals of an excitatory neuron was 1.5 mm, and of an inhibitory neuron was 0.5 mm. The probability of synaptic connection between two nearby excitatory neurons was 0.09. Each excitatory neuron sent a straight 12 mm long myelinated axon to a randomly chosen distant part of the sphere, and its collaterals spanned an area of radius 0.5 mm, as illustrated in Figure 1. Each excitatory neuron innervated 75 local postsynaptic targets chosen randomly from those within a circle of radius 1.5 mm, as well as 25 distant postsynaptic targets chosen from those within a circle of radius 0.5 mm. Each inhibitory neuron innervated 25 randomly chosen neurons within a circle of radius 0.5 mm.

#### Spike Propagation Velocity

The conduction velocity for the myelinated axons was 1 m/s, in accord with the value experimentally measured for myelinated cortico-cortical fibers in awake adult rabbits (Swadlow, 1994). In contrast, the conduction velocity of non-myelinated local axonal collaterals was only 0.15 m/s, in accord with values reviewed by Waxman and Bennett (1972); corresponding delays could be as long as 10 ms.

### Neuronal Dynamics

Each model neuron is described by the following equations (Izhikevich, 2003):

The variable *v* denotes the membrane potential of the neuron and *u* denotes its recovery variable, e.g. the activation of K current. The part 0.04*v*^{2} + 5*v* + 140 was obtained by fitting the spike initiation dynamics of cortical neurons so that units of *v* correspond to mV and units of time correspond to ms. The variable *I*_{syn} denotes the total synaptic current as explained below.

An illustration of the dimensionless parameters *a*, *b*, *c* and *d* is presented in the top of Figure 2, left. The parameter *a* determines the rate of recovery, and parameter *b* determines the sensitivity of the recovery variable *u* to the membrane potential *v*. The parameter *c* determines the after-spike reset value, which depends on fast high-threshold conductances. Similarly, the parameter *d* determines the after-spike reset of the recovery variable, which depends on slow high-threshold conductances.

Values of the parameters corresponding to different firing patterns are presented in the top of Figure 2, right. The lower part of Figure 2 shows membrane potential responses (variable *v*(*t*)) when pulses of DC current (variable *I*) are applied. Despite the simplicity of the model, it reproduces a number of different firing patterns that resemble those observed in the mammalian cortex (Connors and Gutnick, 1990; Gibson *et al.*, 1999).

In the model each excitatory neuron has (*a*, *b*) = (0.02, 0.2). To achieve heterogeneity (so that different neurons have different dynamics), the excitatory neurons were assigned (*c*, *d*) = (–65, 8) + (15, –6)*r*^{2}, where *r* is a random variable uniformly distributed on the interval [0, 1]. Thus, *r* = 0 corresponded to a regular spiking (RS) cell, and *r* =1 corresponded to a chattering (CH) cell. We use *r*^{2} to bias the distribution toward RS cells. Each inhibitory cell had (*a*, *b*) = (0.02, 0.25) + (0.08, –0.05)*r* and (*c*, *d*) = (–65, 2). This heterogeneous distribution is not intended to model any specific cortical area, but rather it reflects a generic principle that the majority of pyramidal neurons in mammalian neocortex are of the RS type (Steriade *et al.*, 2001). We have explored other distributions and found that they did not significantly affect the results.

To avoid numerical instabilities associated with fast spiking activity, each neuron was simulated with a time step of 0.5 ms using the first-order Euler method. The time step used for simulation of synaptic dynamics was 1 ms. The model was coded in C++ programming language. On average, it took 60 s for a 1 GHz Pentium desktop PC to model 1 s of simulated neuronal activity (model time).

### Synaptic Dynamics

#### Input

Each neuron received cortico-cortical synaptic input *I*_{syn} described below and a superthreshold noisy input generated by a Poisson point process with mean firing rate 0 < *q* < 2 Hz, where *q* is a parameter. Most of the simulations were carried out with *q* = 1. This input is the only source of noise in the model.

#### Short-term Depression and Facilitation

We use the phenomenological model by Markram *et al.* (1998) to simulate short-term synaptic plasticity:

where *R* is a ‘depression’ variable, *w* is a ‘facilitation’ variable, *U*, *F* and *D* are parameters that have been measured for a number of cortical neurons (Markram *et al.*, 1998; Gupta *et al.*, 2000), and δ is the Dirac delta function. The product *R*(*t*)*w*(*t*) is the fractional amount of neurotransmitter available at time *t*. Each firing of a presynaptic neuron, occurring at time *t** _{n}*, decreases the ‘depression’ variable R by Rw, and increases the ‘facilitation’ variable

*w*by

*U*(1 –

*w*). Excitatory neurons have

*U*= 0.5,

*F*= 1000 and D = 800. These values correspond to synapses exhibiting depression. Inhibitory synapses have

*U*= 0.2,

*F*= 20 and

*D*= 700. This corresponds to depression synapses classified as type F2 by Gupta

*et al.*(2000).

#### Synaptic Conductances

The total synaptic current of the *i*th neuron is

where *v* is its membrane potential, and the subscript indicates the receptor type. Each conductance *g* (here we omit the subscript for the sake of clarity) has first-order linear kinetics:

with τ = 5, 150, 6 and 150 ms for the simulated AMPA, NMDA, GABA_{A} and BABA_{B} receptors, respectively (Dayan and Abbott, 2001).

The ratio of NMDA to AMPA receptors was set to be uniform at a value of 1 for all excitatory neurons. Thus, each firing of an excitatory presynaptic neuron *j* increases *g*_{AMPA} and *g*_{NMDA} by *c*_{ij}*R*_{j}*w** _{j}*, where

*R*

*and*

_{j}*w*

*are the short-term depression and facilitation variables as above, and*

_{j}*c*

*is the weight of synaptic connection from neuron*

_{ij}*j*to neuron

*i*. The ratio of GABA

_{B}to GABA

_{A}receptors is taken to be 1 for all inhibitory neurons, so each firing of an inhibitory presynaptic neuron increases and by

*R*

_{j}*w*

*.*

_{j}#### Long-term Synaptic Plasticity

In the initial state of the system, all excitatory synaptic weights have the same value, *c** _{ij}* = 0.1. Synaptic dynamics consisted of passive, activity-independent changes and active spike-timing-dependent changes.

The dynamics of passive change of the synaptic weight *c** _{ij}* from neuron

*j*to neuron

*i*are described by the second-order linear equation

where *a* describes slow, activity-independent increase of synaptic weight. The value *a* = 10^{–6} results in augmentation of the strength of all synapses even when the postsynaptic neuron *i* is quiescent. The small parameter 10^{–4} causes the synaptic plasticity in the model to be slow and long-term.

Each synaptic weight, *c** _{ij}*, is updated according to the STDP rule depicted in Figure 3. Depending upon the order of firing of pre- and postsynaptic neurons, we increase or decrease the derivative, . Instantaneous changes of induce slower changes of the synaptic weight

*c*

*according to the equation above. Thus, the STDP rule in our implementation affects the rate of change of synaptic weight rather than the weight itself. Positive values of result in slow potentiation; negative values of result in slow depression. We keep synaptic weights in the range [0, 0.5] by manually clipping*

_{ij}*c*

*.*

_{ij} Unlike an implementation of the STDP rule that considers all possible pairs of pre- and postsynaptic spikes (Song *et al.*, 2000), we assume that only the earliest-neighbor pairs, which are connected by vectors in Figure 3, contribute to plasticity. Firing of a presynaptic neuron *j* at time *t* decreases by , where *t** _{i}* is the time of the last spike of the

*i*th postsynaptic neuron. Thus, the last firing of a postsynaptic neuron overrides the effect of previous firings. Similarly, firing of a postsynaptic neuron

*i*increases by , where

*t*

*is the time of the arrival (after axonal conduction delay) of the last spike of the*

_{j}*j*th presynaptic neuron. That is, the last firing of a presynaptic neuron overrides the effect of previous firings. Such an implementation is consistent with the empirical data reported by Sjostrom

*et al.*(2001), and it has many useful properties (van Rossum

*et al.*, 2000; Rubin

*et al.*, 2001; Izhikevich and Desai, 2003).

#### Spike-timing in Neuronal Groups

In Figure 4 we describe an algorithm for identifying neuronal groups.

Step 1: We find an excitatory anchor neuron (A) having two or more strong connections to other local-circuit excitatory neurons (B, C, D). Here ‘strong’ means that the synaptic weight *c** _{ij}* is within 5% of the maximal value 0.5; see Figure 6. Considering only strong connections ensures that firing of the anchor neuron A enhances the probability that neurons B, C and D fire in the order determined by the conduction delays depicted in Figure 4.

Step 2: We determine the common local-circuit postsynaptic targets of B, C and D. If there are none, then we discard A and proceed to another anchor neuron. If, on the other hand, there are common postsynaptic neurons, e.g. C and D, both synapsed on E and F, then we keep A. Then we check whether the axonal conduction delays and synaptic delays are such that the spikes from C and D arrive at the postsynaptic targets simultaneously (±2 ms). We refer to such delays as being convergent or matching delays. In this illustration, they do not converge simultaneously on E, but do so on F. Similarly, the conduction delays from B and C converge on G. Using the information about the strength of EPSPs (see Fig. 6), we determine whether the input to F and G is superthreshold. If it is, then we consider F and G a part of a group.

Step 3 is the same as step 2 except that we consider common postsynaptic targets of B, C, D, F and G. Neurons I, J and L will be added to the group, while H and K will be discarded because of their non-matching delays. By repeating these steps until either there are no more common postsynaptic targets with matching delays, or all input strengths are subthreshold, we can resolve the neuronal group contributed to by the neuron A.

## Results

For convenience, we review the salient features of the model: It consisted of 100 000 neurons exhibiting known types of spiking and bursting behavior as summarized in Figure 2. The model had 8.5 million excitatory and inhibitory synaptic connections exhibiting both short-term depression and facilitation, as well as long-term STDP. The neurons were arranged randomly on the surface of a sphere, in an arrangement that preserved some of the important ratios found in the mammalian cortex. The ratio of excitatory to inhibitory cells was 4/1, and the probability of a local-circuit connection between any two nearby excitatory neurons was 0.09. Local and global cortico-cortical reentrant connections were mediated by axons having finite conduction velocities giving rise to significant conduction delays ranging from 1 to 12 ms. As described above, these delays played an important role in determining the fine temporal structure of neuronal firing.

### Collective Behavior

Initially, all synaptic weights had equal non-zero values, and the neurons received noisy superthreshold input. The model was allowed to run for 1 day (model time) to allow weights to arrive at asymptotic values. The final state achieved after this one-day simulation was used as the initial state for all subsequent experiments, except the one in Figure 10.

#### Rhythms and asynchronous firing

Figure 5*A* illustrates how the collective dynamics of the model depended on the magnitude *q* of the superthreshold noisy input (see Materials and Methods). The averaged firing rate of all excitatory neurons and the magnitude of the fluctuation in this rate as a function of the parameter *q* are both plotted. When *q* was small (infrequent input), the rate was small and it fluctuated significantly, indicating that the network dynamics exhibited low-frequency oscillation, as illustrated in Figure 5*B*: the activity alternated between quiescence and synchronous bursts. Such rhythmic synchronized behavior is ubiquitous in systems consisting of relatively small numbers of neurons, such as cortical slabs (Timofeev *et al.*, 2000; Bazhenov *et al.*, 2002) or cultures (Latham *et al.*, 2000a,b; Keefer *et al.*, 2001). Increasing *q* resulted in an increase of the firing rate and a decrease of the magnitude of its fluctuation. Values of *q* ≥ 1 Hz (frequent input) resulted in high-frequency, asynchronous, nearly Poissonian neuronal firings, illustrated in Figure 5*C*. In the rest of this paper we use *q* = 1 Hz, which gives rise to inherently stable network dynamics. On average, excitatory neurons fired 10 spikes/s; one spike could be attributed to the noisy input; the remaining spikes were due to synaptic input from the other cortical neurons. Although the number of inhibitory neurons was smaller than that of excitatory neurons, their firing rate was proportionally higher, resulting in an approximate balance of excitation and inhibition in the network.

#### Synaptic Weights

In Figure 6, left, we show the steady-state distribution of excitatory synaptic weights *c** _{ij}*. The bimodal distribution arises from our implementation of the STDP learning rule (Song

*et al.*, 2000) and the fact that the weights are cut off at a maximal value of 0.5. In Figure 6, right, we show the distribution of the corresponding EPSPs. Due to short-term depression and facilitation and the heterogeneity of intrinsic neuronal properties, this distribution is not bimodal, but rather is skewed toward smaller values. The majority of synapses are weak (median EPSP amplitude < 0.1 mV), and only a small percentage produces EPSPs > 5 mV. Weak synapses play an important role in providing background synaptic drive to neurons, and they control the global state of the network, e.g. its rhythmicity. Strong synapses determine the fine temporal structure of neuronal firing, and they are responsible for the spike-timing dynamics leading to the formation of neuronal groups.

### Spike-timing in Neuronal Groups

In Figure 7*A* we show the spiking activity of a small subset of pyramidal neurons in the model. As one can see in the magnifications in Figure 7*B*, there are certain firing patterns that appeared repeatedly in the spike train. Within each pattern, neurons did not fire simultaneously. Instead spikes occurred predictably after fixed temporal intervals, i.e. they were time-locked to each other. Each firing pattern corresponds to an activation of a neuronal group — a set of dynamically interactive neurons having strong synaptic connections with matching conduction delays. The two distinct firing patterns in Figure 7*B* correspond to activation of the ‘blue’ and ‘red’ neuronal groups whose spatial layout is depicted in Figure 7*C*. The causal relationship between spikes and the convergence of delays is plotted in the space-time coordinate system in Figure 7*D*.

Using statistical analysis of spike trains (Oram *et al.*, 1999), we identified many such neuronal groups, although this daunting computational task could not be extended because of the combinatorial explosion. An easier way to discern the formation of spike-timing neuronal groups in the model is to use the available information about the synaptic strengths and conduction delays, which we describe next.

#### Identification and Characterization of Neuronal Groups

Because we know the synaptic weights and conduction delays of all connections within the simulated network, we can determine the order of firing of presynaptic neurons that makes any given neuron fire. Consider, for example, two neurons, X and Y, that have strong synaptic connections to neuron Z, but with different conduction delays, say 10 and 4 ms, respectively. If X fires 6 ms before Y fires, their spikes would arrive at Z simultaneously, thereby evoking a potent postsynaptic response. In contrast, if X and Y fire at a different time interval, e.g. synchronously, then their spikes arrive at Z at different times, and they rarely cause Z to fire. Firing at this different time interval could, of course, cause some other neuron in the network to respond. By considering pre- and postsynaptic neurons one-by-one, this method allows us to determine the structure and spike-timing in each neuronal group in the simulation (see Materials and Methods).

Applying this procedure, we assigned a separate color to each neuronal group and plotted the groups on the cortical surface in Figure 8. Neuronal groups in the model showed quite complicated geometry and firing patterns. Histograms in Figure 9 show that most groups contained fewer than 100 neurons, and they fired a spatiotemporal pattern spanning 25 ms or longer.

#### Growth

Initially, the excitatory synapses have equal strengths (*c** _{ij}* = 0.1 for all

*i*and

*j*) but are connected by different conduction delays that depend on the distance between neurons. Depending upon the pattern of conduction delays, STDP results in depression of some synapses and potentiation of others. As a result of such plasticity, neurons self-organize into groups. The spontaneous appearance of neuronal groups occurs quite quickly, as shown in Figure 10, for the interval 0 to 15 min of model time. After initial growth, the number of groups levels off and subsequently fluctuates within a narrow range, suggesting that a dynamic equilibrium is reached in which new groups appear at the same rate that old groups are lost.

We randomly shuffled all synaptic weights at model time 1000 min. One such experiment is shown in the right-hand side of Figure 10. The number of groups dropped from approximately 900 to 80. Even these 80 groups disappeared within a minute, and after another minute new groups appeared. Thus, only a small number of groups would be generated by chance in a randomly connected network of 100 000 neurons, and these random groups would not persist. This reinforces the conclusion that the majority of the groups in the model were activity dependent and needed STDP to appear and persist.

#### Lifespans

The persistence of each neuronal group in the model was monitored every second during a 24 h simulation in model time (more than a month in real time). In Figure 11 (gray bins) we depict the distribution of the lifespans of neuronal groups. Although many groups had a short lifespan — less than a few minutes — a significant number (almost 25%) of groups persisted during the entire duration of the 24 h experiment.

#### Persistence in the Presence of Synaptic Turnover

There is growing evidence that dendritic spines and synapses in adult brain can be rapidly formed and eliminated with a mean half-life of the order of 4 months (rat barrel cortex: Trachtenberg *et al.*, 2002) or 12.5 months (mouse primary visual cortex: Grutzendler *et al.*, 2002). To investigate the robustness of neuronal groups in the simulation in the face of synaptic turnover, we perturbed all synaptic weights one-by-one: every 10 ms a randomly chosen excitatory synapse, which had not been perturbed before, was reset to zero ( ). Since there are 8 million such synapses in the model, the synaptic turnover time was <24 h, i.e. within this period of time, all synaptic weights had been reset to zero at least once. The distribution of the lifespans of neuronal groups during such an experiment is plotted in Figure 11 (yellow bins). Comparing yellow bins (synaptic turnover) and gray bins (no turnover), one sees that there is essentially no difference between them. Moreover, both bins consist of essentially the same groups, implying that these groups are quite robust.

To account for this robustness it is important to note that synapses were reset rather gradually during the 24 h period. Synaptic connections within a neuronal group consisting of 100 neurons were perturbed on average only once every 15 min, providing sufficient time for activity-dependent self-repair and recovery.

### Effects of Patterned Input

It has been found that persistent alterations in patterns of sensory input to vertebrate cortices can lead to alterations in cortical receptive field properties (Merzenich *et al.*, 1983). We used the model to study how coherent sensory stimulation affects neuronal groups, and the formation of receptive and projective fields. Here, the receptive field of a neuron is the set of stimuli that make the neuron fire. Conversely, the projective field of a stimulus is the set of neurons that fire in response (Lehky and Sejnowski, 1988). We choose two circles of radius 0.5 mm (black circles on the surface of the spheres in Fig. 12) and deliver superthreshold stimuli to all neurons within that circles every 500 ms. Thus, twice a second, the neurons within the circles fire synchronous spikes.

Initially, the projective fields are small. In Figure 12*A*, we use two colors to plot all neuronal groups that have the anchor neuron inside the stimulation circles. Whenever the anchor neurons are stimulated, the groups are activated, and neurons in the groups span the projective field. Thus, all marked neurons have one of the circles in their receptive fields. Neurons that have both circles in their receptive fields are depicted as blue dots.

In Figure 12*B* we show the results of stimulation of the left circle. Coherent firings of neurons within and around that circle propagate to neighboring neurons, thereby recruiting new cells into the existing groups or creating new neuronal groups. As a result, the size of the red projective field grows. After only 5 min (model time) of continuous 2 Hz superthreshold stimulation, the number of neuronal groups originating from the circle is increased by a factor of 60. That is, more neurons are in the projective field of the stimulated circle. Interestingly, the size of the green projective field and the total number of all neuronal groups on the sphere doubles. Nevertheless, the network activity shows a remarkable stability — neurons continue to fire with the same mean firing rate (<10 Hz).

In Figure 12*C* we show the result of synchronous stimulation of both circles. The size of the green projective field increases by a factor of 100, and there are many neurons, marked by blue dots, that have both circles in their receptive fields. The existence of these shared neurons indicates that red groups and green groups cooperate when they are activated simultaneously. Finally, in Figure 12*D* we stimulate both circles asynchronously. Even though the green projective field is still greater than in Figure 12*A* by the factor of 10, it is much smaller than in Figure 12*C*, indicating the presence of a competition between red and green neuronal groups with predominance of the red group. The number of shared neurons is only 10% of that in Figure 12*C*, also supporting the conclusion that asynchronously activated groups compete.

## Discussion

In this paper, we describe simulations of a detailed neuronal network constructed in order to study the collective behavior and self-organization of spiking neurons to form neuronal groups (Edelman, 1987, 1993). Although the model lacks cortical layers (and hence differs from some of our previous detailed models, e.g. Lumer *et al.*, 1997a,b), it preserves many properties characteristic of a mammalian cortex. In particular, neurons in the model exhibit a variety of known types of spiking and bursting patterns; These neurons have receptors with AMPA, NMDA, GABA_{A} and GABA_{B} kinetics, and the synaptic connections among them exhibit both short-term and long-term plasticity. The most critical features governing the dynamics of the model are spiking neurons, axonal conduction delays and STDP. To the best of our knowledge, these three features have not previously been incorporated together in a large-scale mathematical model constructed to shed light on cortical properties.

### Spike-timing-dependent Plasticity

STDP is a novel form of Hebbian synaptic plasticity that has received much experimental and theoretical attention in recent years (Markram *et al.*, 1997; Bi and Poo, 1998; Debanne *et al.*, 1998; Feldman, 2000; van Rossum *et al.*, 2000; Sjostrom *et al.*, 2001; Froemke and Dan, 2002F; Izhikevich and Desai, 2003). According to STDP, the temporal order of presynaptic and postsynaptic spikes determines whether a synapse is potentiated or depressed (Fig. 3). If a presynaptic neuron fires before a postsynaptic neuron, the synapse is potentiated. The smaller the lag between the spikes, the stronger the potentiation. Conversely, the reverse order of firing results in depression of the synapse. We assume that all synapses in the network undergo STDP all the time. The question addressed in the present study is how this plasticity might lead to self-organization in a model in which synapses are considered as populations (Edelman, 1987, 1993).

### Delays

There is a causally important interplay between spike timing and conduction delays, which can be tens of milliseconds in mammalian neocortex (Swadlow, 1994). If spikes arrive at the postsynaptic cell at drastically different times, then the simultaneous firing of presynaptic neurons is not the best mode to fire a given postsynaptic cell (Bienenstock, 1995; Miller, 1996a,b). In order to maximize the probability of postsynaptic response, presynaptic neurons should fire in certain temporal patterns determined mainly by conduction delays (illustrated in Fig. 7*D*). The same neurons firing with different patterns may elicit responses in different postsynaptic cells.

The interplay between delays and STDP gives rise to the spontaneous formation of neuronal groups, i.e. small collectives of neurons having strong connections with matching conduction delays, as in Figure 7. and exhibiting time-locked but not necessarily synchronous spiking activity.

### Neuronal Groups

The neuronal groups considered here can be characterized either by their precise firing sequence or by their synaptic organization. That is, each group, if activated, gives rise to a statistically significant repeatable pattern in the raster of spiking data (e.g. Fig. 7*A*), and each group is linked by a pattern of strong synaptic connections with matching or convergent delays, as in Figure 7*D*.

The spontaneous formation of groups is a consequence of the effect of STDP on the distribution of synaptic delays within the network, since STDP tends to potentiate connections having matching delays that contribute to group formation, and to depress those connections that impede or interfere with this process. At times when a group is not active, the neurons comprising it are not silent (see Fig. 7*A*). Instead, they fire relatively uncorrelated Poisson spike trains that contribute neither strengthening nor weakening of synaptic connections among them. Each group must be activated from time to time, however, to reinforce its existing connections. Prolonged nonactivation results in slow degradation of synaptic strengths and demise of the group. This demise can occur either because some strong connections are depressed, or because some weak connections are potentiated and interfere with the firing pattern of the group as a whole. This process of group dissolution can be accelerated by incongruent patterns of spiking input from other active groups. In this way, most neuronal groups in the network compete with one another for prolonged survival.

### Self-organization

As is evident from Figure 6, the majority of EPSPs are weak, with a mean and a median <0.1 mV. These synapses control the overall state of the network: its rhythms, mean firing rates, etc. The existence of neuronal groups in the model relies on fewer but stronger synapses that are capable of producing sizable EPSPs.

The number of different ways that strong synaptic connections can be arranged to make neuronal groups through spike timing in even a small set of neurons is astronomical. There are, however, far more ways of arranging strong connections so that the neurons fire incoherently and do not form a spike-timing neuronal group (see ‘synaptic shuffling’ in Fig. 10). That is, if the synaptic connections lack matching delays, firing of the neurons will not produce a repeatable spatiotemporal pattern. Thus, we find it significant that neurons do spontaneously self-organize into groups with coherent firing patterns, given the mechanism of STDP applied here.

There are two major ways according to which a small set of neurons might self-organize into a group. First, certain neurons might receive a persistent coherent input that makes them fire repeatedly with a certain spatiotemporal pattern. Among the astronomical number of all possible synaptic organizations, only a few might be consistent with that firing pattern. STDP would potentiate this synaptic organization, effectively excluding others. For example, the pattern of firing of neurons C, H, F and D in Figure 7*D* potentiates red synaptic connections because the spikes from neurons H, F, and C arrive at neuron D just before it fires. A different pattern of firing of the same neurons may potentiate some other synaptic connections and depress red ones.

Alternatively, random firing of the neurons might result in some synaptic connections being potentiated and thus selected by chance, while others are depressed. This initial potentiation will introduce weak correlations of neuronal firing, produce even greater potentiation, and favor some synaptic organizations over others. Once groups are formed, however, their persistent activation may provide a coherent spatiotemporal pattern to local and distant postsynaptic targets and promote formation of neuronal groups at these locations. Our simulation confirms this. Figure 10 shows that there is a rapid increase in the number of neuronal groups with time. Figure 8 illustrates that the groups cover the entire ‘cortical’ surface.

We find that, once formed, many groups (∼25%) are quite stable and can persist for many hours of simulation time (see bin ‘24 hours’ in Fig. 11). Even synaptic turnover, which changed all synapses in a random order during a 24 h period, did not affect the stability of spike-timing neuronal groups (yellow bins in Fig. 11). Repetitive activation of neuronal groups combined with STDP resulted in a self-reconstructing behavior of the model.

### Competition

Many groups are reentrantly connected via either local-circuit connections or global corticocortical connections. Since each group fires a spatiotemporal pattern of spikes that is consistent with its own synaptic organization, activation of one group may distort the firing pattern of other groups, thereby leading to the depression of their internal synaptic connections via STDP.

We can illustrate competitive interactions by referring to the groups labeled red and blue in Figure 7. Because these groups are spatially overlapping, their constituent neurons are anatomically interconnected by numerous intracortical links (not shown in the figure for the sake of clarity). When a neuron in either of the two groups fires, it sends input to neurons in the other group that would affect their firing. This exchange of input acts to disrupt the distribution of synaptic strengths that maintain the firing patterns characteristic of the two groups. If, for example, activity in the blue group in Figure 7 causes neuron B to fire repeatedly before neuron A, this would act to depress the synaptic connection from A to B and lead to dissolution of the red group. In the simulation, the occasional activation of the red and blue groups produced only a weak mutual disruption of intragroup synaptic strengths that allowed both groups to be maintained.

In contrast, persistent coactivation of two or more neuronal groups can have a profound effect, as we illustrate in Figure 13. Activation of either a red or blue neuronal group shown in the figure may not have a significant effect on the yellow group. However, time-locked coactivation of the two groups can project an overwhelming temporal input to the yellow group and make its neurons fire with a pattern that is inconsistent with its own internal synaptic organization. As a result of STDP, synaptic connections are depressed and the yellow neuronal group is destroyed. Simultaneously, coactivation of the red and blue groups might promote the birth of a new neuronal group whose synaptic organization (green in Fig. 13*C*) is consistent with the pattern of spatiotemporal input coming from the two groups.

### Reentry

In this simulation, as in the developing vertebrate brain, neuronal groups do not arise and function independently; they are embedded within a larger network of groups reciprocally linked to one another by long-range excitatory interconnections. It is the ongoing process of massively parallel reciprocal signaling among local neuronal groups, referred to as reentry, that gives rise to the integration of the system as a whole (Edelman, 1987; Tononi *et al.*, 1992).

Since, in the simulation, groups arise independently, their firing patterns are initially uncorrelated with one another and signals exchanged between groups have little effect on intra-group size, neuronal composition or lifetime. Since STDP is a property of synapses mediating reentrant signaling in the simulation, functional connectivity will be enhanced between groups sharing coherent patterns of activity. These patterns may arise initially by chance or they may be the consequence of shared input. The strengths of inter-group connections in the simulation grow to such a degree that reentrant activity can either disrupt or enhance the intra-group firing patterns that determine group structure and longevity. Groups therefore can become competitive or mutually reinforcing, and local group properties are both shaped by and contribute to overall network activity.

In this way, the process of reentry imposes global constraints on the process whereby groups resolve spike-timing conflicts: neuronal groups whose firing patterns conflict with groups to which they are interconnected either adapt their properties or die out, while those whose outputs are consistent with those in other regions are strengthened. Thus reentrant connections establish and maintain cooperation among different regions of the network. This has inportant consequences for information processing by reentrant spiking networks. For example, reentry allows the sensory cortices to unify neuronal signals arising in different sense organs from a single stimulus, in this way resolving the binding problem in a fashion consistent with previous modeling studies (Tononi *et al.*, 1992).

### Neuronal Groups and Synfire Chains

The concept of a neuronal group governed by spike timing as considered here is quite different from the notion of a synfire chain proposed by Abeles (1991, 2002), Diesmann *et al.* (1999) and Aviel *et al.* (2003), although both concepts involve temporally connected spiking. Synfire chains consist of pools of neurons connected sequentially so that synchronous firing of all neurons in one pool propagates to the next pool without significant temporal distortions. Having heterogeneous axonal conduction delays would be devastating for a synfire chain, since this would lead to instability and breakdown of synchronous activity. In contrast, neuronal groups do not have any identifiable pools, but actually rely heavily on the heterogeneity of conduction delays as in synfire braids (Bienenstock, 1995). This heterogeneity, which is consistent with the TNGS (Edelman, 1987, 1993), tends to reduce local synchronous firings and makes local spiking dynamics time locked and asynchronous.

### Statistics of Spike Rasters

It is relatively easy to find neuronal groups using anatomical data from the model and the algorithm presented in the Materials and Methods section. Once the firing pattern of any group is determined, one can use the pattern as a template and scan the spike train to determine when the group is activated, as we do in Figure 7*A*. This method is free from the argument (Oram *et al.*, 1999; Baker and Lemon, 2000) that precise firing sequences result from random spike trains with co-varying firing rates.

In contrast, finding groups using only statistical analyses of spike trains is practically impossible even in the model, let alone *in vivo*. Indeed, only 250 out of 1000 spike-timing neuronal groups live long enough to reveal themselves in the spike train (see Fig. 11). Since on average in the model there are fewer than 30 neurons per such a group (see Fig.9), fewer than 7500 neurons (<7.5% of all the neurons in the model) are part of any such group. A simple statistical calculation suggests that one would need to record from 866 neurons simultaneously to achieve a probability of 0.5 of identifying three neurons belonging to the same neuronal group in the network. Thus, spike rasters of thousands of neurons might be needed to reveal even the smallest spike-timing group in a functioning neocortex. Such a large-scale recording over a long time period is infeasible at present.

## Conclusion

We have simulated a network of neurons having detailed spiking and synaptic dynamics, STDP and realistic axonal conduction delays. Despite the lack of coherent input in the model, neurons self-organized into groups and repeatedly generated patterns of activity with a millisecond spike-timing precision. The properties of the model allowed us to explore the spike-timing dynamics of neuronal groups, in particular their formation, activity-dependent persistence, competition, cooperation and responses of receptive fields to coherent input. We hypothesize that the appearance of such neuronal groups is a general property of reentrantly connected spiking networks with axonal conduction delays and STDP.

Niraj S. Desai, Anil K. Seth, John R. Iversen, Elisabeth C. Walcott, Yuval Aviel, George Reeke, and Jeffrey Krichmar read early drafts of the manuscript and made a number of valuable suggestions. This research was supported by the Neurosciences Research Foundation.