David E. Huber (dhuber@indiana.edu)

Psychology Department

Indiana University

Bloomington, IN 47405

A long standing problem in theories of information processing is the lack of adequate mechanisms for binding together separate components of a stimulus or memory. Traditionally, theorists have proposed high-level abstract models to deal with this "binding problem" without appealing to the neural substrate of processing (Biederman, 1987; Carpenter & Grossberg, 1988; Treisman & Gelade, 1980). Recent evidence demonstrates that the temporal correlation of neural firing patterns (synchrony) exists in many species under conditions of normal behavior and is stimulus specific in its elicitation (for a review article see Singer & Gray, 1995). Accordingly, some theorists propose that synchrony is the neural mechanism underlying binding within their models (Grossberg & Somers, 1991; Hummel & Biederman, 1992; Treisman, 1996).

In support of these ideas, numerous simulations and proofs demonstrate the ease with which coupled oscillators readily achieve a synchronous state. Synchrony has been demonstrated for leaky integrate and fire models of neurons which couple by discrete activation changes (Gomez & Budelli, 1996; Mirollo & Strogatz, 1990), simple sinusoidal oscillators coupling through phase pulling (Lumer & Huberman, 1992), two variable relaxation oscillators with continual coupling (Grossberg & Somers, 1991), and models implementing multiple ion currents through Hodgkin-Huxley style gating terms (Golomb, Wang, & Rinzel, 1994; Demir, Butera, DeFranceschi, Clark, & Byrne, 1997).

The wide range of models and conditions under which synchrony results suggests that synchronous oscillations might be unavoidable. Virtually universal synchrony would not serve a useful function and might predict almost universal epilepsy (for classic work on synchrony and epilepsy see Jasper & Kershman, 1941). In this article I provide evidence that the buildup and propagation of action potentials along nerve fibers substantially limit the range of conditions under which synchrony develops.

Somers and Kopell (1993; 1995) provide mathematical proofs that the group of models termed relaxation oscillators, referring to systems operating with variables on different time scales, synchronize more readily than other models. Their theorems are specific to two variable oscillators such as the Morris-Lecar (Morris & Lecar, 1981) or Fitzhugh-Nagumo (Fitzhugh, 1961; Nagumo, Arimoto, & Yoshizawa, 1962) models, but are relevant to the four variable Hodgkin-Huxley equations (Hodgkin & Huxley, 1952) which the two variable models approximate. Essentially these relaxation oscillators easily synchronize since small impulses cause multiple cells to immediately and synchronously "fire" due to rapid threshold modulations. This situation can arise even between cells of differing frequencies since the fast variable (i.e. Na+ gating) is unchanged with frequency.

In order to demonstrate that inclusion of spatial propagation provides an important constraint on synchrony development, I select one of these relaxation oscillators. Rather than working with the four variable Hodgkin-Huxley equations, I choose the two variable Fitzhugh-Nagumo model. Limiting the situation to two variables allows for a qualitative accounting of behavior through phase portrait analyses.

The Fitzhugh-Nagumo model was independently derived by Fitzhugh (Fitzhugh, 1961) and Nagumo (Nagumo, Arimoto, & Yoshizawa, 1962) from the Hodgkin-Huxley equations by assuming that Na+ gating is instantaneous and lumping K+ gating, leakage currents, and ATP pump action into a single recovery variable. As a simplification of ion currents, it is inadequate for quantitative modeling yet sufficiently captures the dynamics for the present situation.

The assumption for Na+ gating results in a cubic expression for thresholding on membrane potential. This is the fast variable, v, whose actions are dictated by the partial differential equation:

(1)

(2)

One way to analyze the behavior of dynamical systems is with a phase portrait. This is a graph (see Figure 1) representing the change in each variable by a vector (arrow) as a function of the current values for each variable. In this graph there is no momentum and subsequent values for the variables are completely determined by their present values. An aid for interpreting phase portraits is the portrayal of isoclines. These are lines along which one of the variables does not change. Isoclines are derived by setting the partial differential equations equal to zero. For equations 1 and 2, this produces:

(3)

for the membrane potential (v) isocline and:

(4)

for the recovery (w) isocline.

Figure 1 shows the phase portrait when the model is driven with an input of I=.112. This corresponds to a real neuron which is fully activated by synaptic long lasting input. For a range of constant inputs (see Figure 2), the model will display this limit cycle behavior. Limit cycle behavior is characterized by a repelling fixed point. A fixed point is any location where the isoclines cross resulting in no change in either variable. With the exception of starting at the fixed point, a cell placed at any other combination of membrane potential and recovery will ultimately relax onto the gray line representing a path of continual oscillations. This is the limit cycle.

For somewhat lesser or greater values of constant input, the membrane potential isocline (eq. 3) will correspondingly shift downwards or upwards and the cell will enter an attractor state. In an attractor state, a cell placed at any combination of membrane potential and recovery will ultimately end up at the fixed point and oscillations will cease (see Figure 2). A cell at rest (I=0) would be an example of this. The cell will remain at rest indefinitely. However, the cubic thresholding allows that a small temporary input to the cell will cause it to cycle a single time and then return to its resting point. This is in keeping with real neurons at rest with some spontaneous firing due to random fluctuations in input.

In order to measure synchronous behavior or even record the frequency of a cell, an explicit firing criterion is needed. I have chosen the membrane potential corresponding to the right hand peak of the membrane potential isocline. For a cell to reach the right hand branch of the membrane potential isocline, it is necessary to cross this criterion. This is analogous to recording voltage fluctuations in a real neuron and using a voltage criterion for determining the exact moment when a cell fires.

All of the previously mentioned models (including the Fitzhugh-Nagumo model) are expressions of inward and outward currents (i.e. radial currents) for any arbitrary point along the axon. The models disregard current in the direction of the axon (i.e. axial current). Radial currents are due to voltage gated or leakage ion channels and the actions of ATP pumps. In order to relate these radial currents to axial current, a simple two step derivation has been developed by applying conservation of charge and Ohm’s law (Fitzhugh, 1962; Jurisic, 1987). This leads to the following expression:

(5)

where the variable x is position along the axon, c is membrane capacitance (which can by a function of the membrane potential), r is axon radius, R is resistance of the axoplasm, and J_r is radial current. This is a very general expression and any of the previously mentioned models for radial current could be used to replace J_r.

In deriving an expression for radial current, Hodgkin and Huxley (1952) used empirical data obtained with a voltage clamped neuron. With this technique the voltage at each position along the axon is kept at a fixed value through a feedback loop. For this reason the second derivative of v with respect to x is assumed to be zero. Based upon these condition, the Hodgkin-Huxley model and others similar to it are only appropriate for explaining the nature of voltage clamped data. In order to capture the true buildup and propagation of an action potential it is necessary to consider the spatial variable as well.

The equations are kept relatively simple by assuming capacitance is minimally dependent upon membrane potential and is set to one. Similarly, the expression r/2R is set to one yielding the following equation for change in the membrane potential within the Fitzhugh-Nagumo framework:

(6)

The axon hillock is placed at position x=0. This is the only position receiving external input. Every other position remains at rest except for changes initiated by axial currents. The equation for the partial differential of recovery with respect to time is the same as before (eq 2) and is applied separately at each position along the axon. Insulated boundary conditions are assumed and the partial differential with respect to x is solved through a spatially centered difference scheme with step size 1.0. Partial differential equations with respect to time are calculated using a forward difference scheme with fixed step size .25.

The reason for analyzing behavior in terms of period instead of frequency is to compare refractory period to baseline period. Refractory period is assessed by driving a cell with a sinusoidally varying input of the same amplitude as the constant baseline input. The period of oscillation for the sinusoidal input was systematically varied from zero up to a value corresponding the baseline period. The shortest period for constant firing was recorded as a rough estimate of refractory period. For simulations with propagation, firing is determined at the end of the axon (x=19). An oscillation occurring at x=0 is irrelevant unless it is capable of producing an action potential which traverses the entire length of the axon.

For assessing synchrony, a grid of 11 X 11 cells was simulated with different degrees of connectivity. At one extreme, each cell communicated with its 8 immediate neighbors (local connectivity). At the other extreme every cell communicated with every other cell (full connectivity). For local connectivity the ends of the grid were connected to form a torus thus avoiding edge effects. All the cells were driven with a constant baseline input (I=.112 without propagation and I=.25 with propagation). For every simulation, cells were placed at randomly determined positions along their limit cycles to place the grid in a completely asynchronous state. The only means for synchronization was through impulses delivered between cells with each passing of their firing criteria.

As with communication between real neurons, these impulses were exponentially decaying inputs. The non-committal term impulse is used since the effect of one cell on another might be due to neurotransmitter gated ion channels or direct electrical coupling (i.e. gap junctions). For each cell the totality of its input {I} was the baseline input plus the sum of the incoming impulses. For each simulation, these impulses were of a fixed amplitude and exponential decay rate. Individual amplitudes were determined by dividing total amplitude by the number of connected cells. In other words, amplitude was divided by 8 for local connectivity and 120 for full connectivity. This equates the connectivity conditions for a fully synchronized grid. 100 simulations were performed with variations of total impulse amplitude from .1 to 1.0 in increments of .1 combined with variations in average impulse duration from 25 to 250 in increments of 25. Each of the bars shown in Figure 4 is averaged across these 100 amplitude/duration combinations.

Synchrony was measured in the same manner found in the single cell recording literature. Auto-correlograms were computed with bins of 5 time units for every cell and then averaged. Cross-correlograms between every cell and every other cell were calculated and then averaged. After an initial 1000 time unit period, cell firing times were recorded for 2560 time units. Recording two simulations for each condition allowed for the calculation of shift predictors to normalize correlograms for random matching (see Singer & Gray, 1995 for a discussion of all these measures). Synchrony was determined by taking the sum of the zero centered 3 bins in the averaged cross-correlogram and dividing by the zero centered bin of the auto-correlogram. If all the cells fired within 15 time units of one another, this would be measured as perfect synchrony (ratio of 1). It should be noted that this is a measure of global synchrony for the entire grid and would under represent a situation of separate synchronous groups of neurons placed at evenly spaced phase relationships.

For simulations with spatial propagation, every cell had an axon of length 20. Perhaps a more realistic assumption would be to vary axon length with the separate distances for each connection. This would serve to further diminish what turned out to be greatly reduced synchronization.

The four sets of results seen in Figure 4 stem from three conditions without the inclusion of propagation and a fourth with propagation. The immediate condition is the standard Fitzhugh-Nagumo model with immediate delivery of impulses between cells. This condition is similar to much of the previous work on synchrony development. In the delay condition, the basic model is augmented with a delay of 60 time units which is approximately the time it takes for propagation down an axon of length 20. The delay+ref condition is a further augmentation through the inclusion of refractory periods as dictated by the results in Figure 3 for a cell with propagation.

Figure 4. Synchrony results. Each bar is averaged across 100 combinations of impulse amplitude and impulse duration.

Approximating the full scale propagate simulation with the addition of delays and refractory periods induced the same ordering of conditions, yet synchronization is substantially less for the propagate simulation. Explicit delays and refractory periods are inadequate for capturing some of the non-linearities involved in propagation. For instance, refractory period is more appropriately a dynamic construct. Recent failed attempts at propagation can block subsequent action potentials almost as strongly as recent successful action potentials. In some sense there is a dynamic memory to the axon. Another example of dynamic non-linearity is found with delay time. While a continually firing cell has a delay time of 60 time units, a fully recovered resting cell is capable of propagation in 40 time units.

These and other anecdotal accounts are provided by watching two real time computer graphics programs developed to analyze the cells within a phase portrait and a x vs. v/w plot.

Replicating previous work, the grid readily synchronized for immediate excitatory coupling for both identical and non-identical cells. This was regardless of connectivity. In contrast, any sort of inhibitory coupling greatly diminished synchronization in the immediate condition. In retrospect this is not surprising. If there is any distribution in firing times the message by the lead cells to the followers will be to not fire making coherency of the group difficult. The addition of a delay, whether it be explicit or through propagation, puts excitation and inhibition on equal grounds for identical cells.

For non-identical cells, the imposition of a delay causes inhibition to produce greater synchrony than excitation. For variable frequency cells, the spreading out of the group during recovery is avoided with inhibition. Without inhibition the faster cells would recover more quickly and fire earlier, whereas the dynamics of an inhibiting impulse override and the group is quiescent until released from inhibition.

The overall effect for extent of connectivity is a reduction in synchronization with local connections. This patterns hold true across all the variables and is due to increased input variability with local connections. For the case of full connectivity every cell in the grid experiences exactly the same input at every point in time. This is in contrast to the situation for local connectivity which is equated in terms of average input but not variability. For local connectivity each cell is following a unique dynamics as dictated by its particular collection of recent incoming impulses.

In general there was a reduction in synchronization for the non-identical (variable frequency) simulations. This reduction is more pronounced for the three conditions including a delay (delay, del+ref and propagate) with the notable exception of fully connected inhibitory coupling. The previously mentioned interaction between frequency variability and inhibition/excitation accounts for this exception. Inhibition provides a more consistent overriding dynamics for full connectivity since it is guaranteed that every cell will be released from inhibition at the same time.

The biologically plausible scenario for coupling through chemical (i.e. neurotransmitter mediated) synapses is the case of non-identical cells in the propagate condition. For these simulations, synchronization was almost nonexistent for local connectivity. Furthermore, fully connected inhibitory coupling produced greater synchronization than excitatory coupling. This result is in keeping with recent reports implicating the inhibitor GABA in the development of synchrony (MacLeod & Laurent, 1996).

The appropriate scenario for electrical synapses (gap junctions) is locally connected excitatory coupling for non-identical cells in the immediate condition. Gap junction impulses are always excitatory and only between neighbors. This kind of synapse results from direct links between adjacent cell bodies making the simulation of propagation irrelevant. Gap junctions are primarily found in low level sensory areas, such as the retina, where cells are densely packed together. In agreement with studies of cells connected through electrical synapses (Neuenschwander & Singer, 1996), synchronization was sizable for this simulation.

These simulations demonstrate that the buildup and propagation of action potentials are an important consideration in analyzing boundary conditions for synchrony development. Propagation is marginally approximated through the inclusion of explicit delay times and refractory periods, yet only a full simulation with a spatial variable and its resulting non-linearities can account for the diminished range of conditions leading to synchronization. Now that a plausible model has been developed which fails to synchronize except in certain situations, the search for architectures resulting in stimulus specific synchronization can be initiated.

___________________