Read, Heather L. and Siegel, Ralph M. (1996) On the origins of aperiodicities in sensory neuron entrainment. Neuroscience vol 75 no 1, pp301-314.

Acknowledgments: We extend our thanks to Dr. Charles Tresser for discussions of oscillatory system dynamics and to Dr. Leon Glass and Dr. James Chrobak for comments on the manuscript. Supported by NIH R01 EY09223-01, ONR N00014-92-1-0334 and a Henry Rutgers Fellowship.

Abstract

Aperiodic entrainment to rhythmic sensory input was obtained with either a single neuron or an excitatory network model, without addition of a stochastic or "noisy" element. The entrainment properties of primary sensory neurons were well captured by the dynamics of the Hodgkin-Huxley ordinary differential equations with a quiescent resting state or threshold for spike output. The frequency-amplitude parameter space was compressed and aperiodic regimes were small in comparison to those of periodically activated pacemaker like neurons. Transitions between phase-locked and aperiodic entrainment patterns were predictable and determined by the equation dynamics; supporting the contention that some aperiodicities observed ${in}$ ${situ}$ arise from the inherent membrane properties of neurons. When the rhythmically activated neuron was embedded in an excitatory network of Hodgkin-Huxley neurons with heterogeneous synaptic delays, aperiodic entrainment patterns were more frequently encountered and these were associated with asynchronous output from the network. Embedding the rhythmically activated neuron in a network with synaptic delays, greatly reduced the range of entrained spike frequencies. Other biological mechanisms of modifying the entrainment properties and promoting aperiodic entrainment are discussed.

Key Words: sensory, neuron, irregular, synchronous, asynchronous, dynamic, entrainment

Introduction

One of the basic organizing principles of coherent activity in the mammalian nervous system is the genesis of rhythms. Neurons are capable of generating cyclic repetitive output, the inherent frequency of which depends on the complement of membrane conductances underlying its spiking behavior. If a neuron is driven with a sufficiently large amplitude rhythmic membrane depolarization it will generate action potentials (spikes) that occur at a fixed time in relation to the stimulus for a limited range of stimulus frequencies. For example, somatosensory fibers can spike with a fixed time delay following the onset of each cycle of membrane depolarization that occurs with rhythmic mechanical vibration of the skin [45]. Spiking is considered "phase-locked" or "periodically entrained" to the stimulus when there is a fixed temporal relation between the cycle (or phase) of input and the cycle of output (spiking) of a neuron. If the stimulus amplitude is sufficiently large, a somatosensory fiber will generate one spike for every stimulus cycle yielding a periodic entrainment pattern of 1:1 (cycles of input: spike output).

Simple 1:1 entrainment patterns are rare in the mammalian central nervous system, the more common observation being a spike pattern that is "phase-related" to the period of rhythmic input, but not periodic itself. Physiologists generally attribute such $jitter$ or aperiodic entrainment to ion channel and synaptic noise; however, a source of this irregularity that is often overlooked is that arising from the entrainment process or the neural system itself.

Previous physiologic and theoretical studies have focused on the entrainment properties of pacemaker like neurons that spontaneously oscillate at rest [1,8,17,19,21,27,32,38,39]. In these biological models, a resting oscillatory state was created by lowering the external calcium concentrations [1,27] or by injecting depolarizing bias current [8,17,19,21]. In either case, periodic input was able to generate aperiodic spike trains in the driven neuron. Though these model systems were simple, they were biological and there was some debate as to the origin of aperiodic entrainment patterns [19,21]. These aperiodic entrainment states exist in simplified mathematical models which do not incorporate noise [1,17,18,19,21,27,31,32]. Furthermore, these aperiodic entrainment patterns were obtained in a parametric and predictable manner similar to that of periodic entrainment. These findings supported the contention that aperiodicities in biological systems are in part a product of the entrainment process.

While the spontaneous oscillatory state is a natural one for cardiac fibers and for pacemaker like neurons, it is not the natural resting state for a broad spectrum of neurons in the mammalian peripheral and central nervous systems [11,33,40,45]. A variety of sensory neurons including: cutaneous afferents, gustatory neurons, warm fibers and bulboreticular neurons have zero or low spontaneous firing rates [4,11,33,40,45]. The dynamics of such neurons have not been realistically captured with computational studies that assume baseline periodicity prior to entrainment. Furthermore, these studies assumed and required the acquisition of continuous intracellular membrane potential data, a technical feat rarely achieved in the natural state of the neuron (i.e., in awake behaving animals).

The present study seeks to explore in detail the entrainment properties of neurons that are silent prior to sensory or afferent input. Although there is an examination of intracellular membrane potential, most of the study is aimed at an easily obtained physiological measure, in particular the interspike interval (ISI), to facilitate comparison with experimental data. Sinusoid current is used to drive otherwise quiescent Hodgkin-Huxley (HH) type model neurons. In addition, the aperiodic entrainment patterns of a sinusoid driven HH neuron embedded in a population of synaptically coupled HH neurons are examined, again to more closely approximate $in$ $situ$ conditions. The synaptic feedback and delays associated with the population of coupled HH neurons serve to dramatically alter the entrainment properties. We find that the range of behaviors of a driven single neuron can be quite broad and could account for much of the temporal patterns or dynamics of primary sensory neurons. The linking together of these neurons into a small population restricts the total range of mean spike frequencies and extends the range of aperiodic entrainment patterns generated with a broad range of stimulus parameters.

Methods

The HH differential equations were used to model the sodium ($g sub NA$), potassium ($g sub K$) and leak ($g sub L$) conductances underlying the action potential [20]. These equations were originally developed to model squid axonal current and they were employed here to provide a simple model of sensory neuron action potentials. The differential equations for the active membrane state variables and activation/inactivation rate functions were as follows: $ {C sub m}dV over dt = {I sub input - g sub Na m sup 3 h ( V - V sub Na ) - g sub k n sup 4 ( V - V sub K ) - g sub L ( V - V sub L )} [1] $ $ dm over dt = alpha sub m [ 1 - m(t) ] - {beta} sub m {m(t)} $ $ dh over dt = alpha sub h [ 1 - h(t) ] - {beta} sub h {h(t)} $ $ dn over dt = alpha sub n [ 1 - n(t) ] - {beta} sub n {n(t)} $ $ alpha sub n = 0.1{{({{V + 10} over 10})} over {exp {({{V + 10} over 10})} - 1}} beta sub n = 0.125 exp {({V over 80})} $ $ alpha sub m = {2.5{({{V + 25} over 25})}} over {{exp {({{V + 25} over 10})} - 1}} beta sub m = 4 exp {({V over 18})} $ $ alpha sub h = 0.07 exp {({V over 20})} beta sub h = 1 over {exp {({{V + 30} over 10})} + 1} $

The variable $I sub input$ is the input current density (${n A} over cm sup 2$) and V is the membrane potential difference (mV). The rising rate constants $alpha sub m$, $alpha sub h$ and $alpha sub n$ account for sodium channel activation, sodium channel inactivation and potassium channel activation, respectively. The decay rate constants for the respective ion channels are $beta sub m$, $beta sub h$ and $beta sub n$. The maximal sodium, potassium and leak conductances are $g sub Na$, $g sub K$ and $g sub L$ ($ mS over cm sup 2 $), respectively; and the respective ionic reversal potentials are $ V sub Na$, $V sub K$ and $V sub L$. The variable $t$ is time (msec) and $C sub m$ is the membrane capacity ($ {mu F} over cm sup 2 $). Temperature, $C sub m$, $g sub Na$, $g sub K$, and $g sub L$ are $6.3 sup o C$, 1.0 ${mu F} over cm sup 2$, 120, 36 and 0.1 $ mS over cm sup 2$, respectively (see table 1 for appropriate values and abbreviations).

Single model or "HH" neurons are given by a complete set of the equations for a unit area of membrane as outlined above. The computational simplicity of these model neurons allows for extensive parameter mapping of spike entrainment. Runge-Kutta fourth order numerical integration routine (time step of 0.05 msec) is employed to integrate these differential equations [36]. Previous studies of these equations have found qualitatively similar results with lower order methods of integration such as the Euler method [21]. However, the fourth order Runge-Kutta method for stiff differential equations is more suitable for quantitative analysis as it introduces less numerical error [19]. Computations were carried out on IBM RS/6000 workstations with 32 bit precision.

The population model consists of 25 coupled HH neurons in a one-dimensional array. Each neuron within the 25 cell array is coupled to its 5 nearest neighbors on each side via a rapid depolarizing membrane current. The synaptic strengths and delays are the same for all neurons. A Heavyside function was used to model the summed synaptic input current for computation simplicity. The summed synaptic input current to a neuron is modeled according to the following function, $ I sub syn sub i {(t)} = sum from {j=1} to {n sub syn} {S sub ij \(bu ({V sub j {(t - {tau sub ij}) - V sub th}})} {H [{V sub j {(t - {tau sub ij}) - V sub th}}]} [2] $ H(x) is the Heavyside function (i.e. $H(x) = 0$ for $x <= 0$; $H(x) = 1$ for $x > 0$). A synaptic potential occurs when $V(t) - V sub {th}$ is greater then zero. The resulting function resembles the $alpha$ function ($I sub {alpha} = texp(-t/{tau}$)) which is commonly employed to model synaptic currents. Fits of the model synaptic function with an $alpha$ function require time constants (${tau}$) in the range of 0.7 to 1.2 msec. The time constant of the model synaptic potential is near that of summed synaptic potentials observed biologically (e.g., mossy fiber synapse [5]). The maximal synaptic potential amplitude increases in proportion to the level of $V sub th$ which was fixed at 50 mV. The synaptic strength constant, $S sub ij$, is fixed at 0.142 $mS over {cm sup 2}$, which yields a subthreshold synaptic potential in the present model. Thus, a certain degree of synaptic convergence must occur to bring the membrane potential to threshold (see figure 6). The conduction delay is the time (t) required for the action potential to travel from the $i sup {th}$ to the $j sup {th}$ neuron and includes the axonal propagation delay and synaptic delay. The conduction delay between adjacent neurons was fixed at 3 msec. This delay increased linearly with displacement, so that the largest delay was 15 msec for input displaced by five neurons from the source. The number of synaptic inputs to a given neuron ($n sub syn$) varied from 5 to 10 depending on the cell position. Thus, neurons displaced 6 or less neurons from the borders had ${n sub syn} = {n sub dis} -1 + 5$, where $n sub dis$ is the displacement from the border; all other neurons received ten synaptic inputs. Neurons were never connected to themselves via recurrent connections. The parameters $S sub ij$, $tau sub {ij}$, and $n sub syn$ were based on extensive neuronal simulation to provide a population single quiescent neurons that would spontaneously oscillate (Siegel and Read, unpublished observation).

The driving current or $I sub input$ for single and population neurons differed in the following manner. In the single neuron model, the applied stimulus current was varied according to a cosine function, $ I sub input = I sub {cos(t)} + I sub {syn(t)} [3] $ $ I sub cos = A cos ({omega} t) $ Where, $f$ = $omega over {2 pi}$. Thus, for the single neuron model $I sub syn$ was zero and $I sub input$ was equal to $I sub cos$. In the population model, only the central neuron was driven by a sinusoidal current and population synaptic input; hence, $I sub input$ is equal to the sum of $I sub cos$ and $I sub syn$ for neuron 13. For all neurons except the central neuron in the population model, the input current was synaptic (i.e., $I sub input$ = $I sub syn$). The amplitude, A (${{nA} over cm sup 2}$), and $omega$ ($rad/{msec}$) were varied systematically in order to map the amplitude-frequency parameter space.

The interspike intervals (ISI's) were computed as the time between the crossings of 50 mV on the upstroke of the action potential. For all figures, the ISI's have been normalized (nISI) to the driving period (T=$1/f$) yielding a dimensionless quantity (i.e., $nISI$ = $ISI(msec)/{T (msec)}$. Normalization of the ISI's allowed for comparison of ISI's generated with different sinusoid frequencies.

Results

The sinusoidally driven single neuron: periodic patterns

The single neuron entrainment to the period of sinusoid input current was analyzed for a large range of sinusoid amplitudes and periods. Single neuron activity was simulated for 50 sec in order to ensure that the initial transients in activity were gone and the output pattern had stabilized. The spike patterns that occurred in the final 10 sec of a 50 sec simulation were examined for periodicity. Spiking patterns were considered phase-locked if the action potential occurred at a regular phase of the stimulus cycle. The means and standard deviations for normalized ISI's (nISI's) occurring in the 10 sec steady-state period were calculated as numerical measures of entrainment.

Spiking of the single HH neuron was temporally locked to the period of sinusoidal current for most stimulus parameters employed; and phase-locked patterns varied as a function of stimulus parameters. With low stimulus intensities, frequent but periodic spike failures occurred and resulted in "higher order" entrainment ratios. A common higher order entrainment pattern consisted of 3 cycles of sinusoid current per 1 spike or a 3:1 ratio (figure 1A). (Throughout this paper the terminology m:n refers to m cycles of sinusoid input to n action potentials). In figure 1 (A, B, and D) the sinusoid stimulus was plotted (-15 mV offset) along with the neuron membrane potential to allow for visual comparison of the stimulus phase and spiking pattern. Patterns of 4:1, 5:1, 6:1, 8:1, 9:1, 10:1 and 11:1 were observed as well (e.g., 11:1 figure 4A); however, these were not as common as the 3:1 pattern of spiking illustrated in figure 1A. There was no deviation from a 3:1 phase-locking pattern in the range of stimulus amplitudes including 1500 to 1515 (${{nA} over cm sup 2}$) with a fixed sinusoid period of 19.04 msec ($omega$=0.31 $msec sup -1$). In the region labeled "3:1" in figure 1C, the mean nISI per stimulus condition was 3.0 and the standard deviation was zero, reflecting the stability of 3:1 phase-locking.

Abrupt transitions from stable higher order to stable lower order phase-locking patterns took place when the sinusoid amplitude was increased. For example, when the stimulus amplitude was increased from 1515 to 1520 (${{nA} over cm sup 2}$), there was an associated transition from 3:1 to 2:1 phase-locking (e.g., figure 1C, $omega$=0.33 $msec sup -1$). Between the regions of 3:1 and 2:1 phase-locking, there was an abrupt increase in the nISI standard deviation. Several small regions of stable phase-locking were found inside this border region. For example, the most stable pattern in the "border region" shown in figure 1C was the 5:2 pattern (see voltage plot of figure 1B). ISI durations of exactly 2 and 3 times the stimulus period were interleaved in the 5:2 pattern. Accordingly, the mean nISI (plotted as a line in figure 1C, arrow) was 2.5 for the 5:2 pattern. A small nISI standard deviation (0.5) was observed with the 5:2 pattern reflecting the stable alternating sequence of ISI's. Thus, the 5:2 pattern was a merging of the stable patterns in the bordering regions of 3:1 and 2:1 phase-locking (figure 1). Similarly, an intermediate pattern of 3:2 was found between 2:1 and 1:1 phase-locking regions of the parameter space. The 3:2 was composed of alternating ISI's of 2 and 1 times the sinusoid period and a mean nISI of 1.5. Thus, the intermediate 3:2 pattern was a merging of stable patterns as was the case for 5:2 patterns. The general trend was for fewer spike failures to occur with higher amplitude stimuli.

The order of transitions between entrainment patterns was conserved for different stimulus periods. For example, the sequence of 3:1, 5:2, 2:1 patterns elicited with an increase in stimulus amplitude was the same for stimulus periods of 29.9, 20.94 and 19.04 msec. However, the range of stimulus amplitudes eliciting a given pattern varied with stimulus period. Thus, the range of stimulus amplitudes eliciting mean nISI's of 3, 2.5 and 2 (phase-locking of 3:1, 5:2 and 2:1, respectively) was smaller with a stimulus period of 29.9 msec (figure 4C) as opposed to 19.04 msec (figure 1C).

The order of patterns observed with different stimulus conditions was predictable according to a simple additive rule. For example, entrainment ratios of 3:1 and 2:1 add to generate a ratio of 5:2. Additivity of neighboring phase-locking patterns was a general property of the sinusoid driven HH neuron. This additive property follows the general rule that regions of $({m sub 1} + {m sub 2}):({n sub 1} + {n sub 2})$ occur between neighboring patterns of ${m sub 1}:{n sub 1}$ and ${m sub 2}:{n sub 2}$ (cycles of input:spikes). Phase-locking ratios including 7:2, 8:2, 8:7, 10:9 and 9:8 were observed sandwiched between appropriate neighbors (e.g., 7:2, figure 4C). This adding rule is not fortuitous; it is described by the "Farey sum" in number theory and can be analytically demonstrated to describe the parametric sequencing of temporal patterns in certain mathematical and physical oscillating systems [10,12,32,39,43]. The sequence of summating patterns are called "Farey sequences".

The rigid order of transitions between phase-locking patterns was conserved throughout the frequency-amplitude parameter space. A high resolution map of the mean nISI's for a large range of parameters was generated to illustrate interactions between the amplitude and frequency of sinusoidal stimuli. Stimulus frequency or $omega$ ($2 pi{f}$) was varied from 0 to 1 ${msec sup -1}$ in steps of 0.005 ${msec sup -1}$. Stimulus amplitude was varied from 2 to 80 in steps of 1 ( ${mu A} over cm sup 2$ ). The steps in amplitude used to map large portions of the parameter space (${mu A} over cm sup 2$; figure 2A) were 2000 times greater than those used to examine Farey sequences in detail (0.5 ${nA} over cm sup 2$; figures 1, 3 and 4). Mean nISI's were assigned to a repeating rainbow color spectrum with mean nISI's of 1, 2 and 3 labeled with fuchsia, yellow and green. The exact values of the color spectrum are labeled on the color bar above the graph in figure 2A. As known for other excitable or threshold systems, the different phase-locking regions all converged towards a common low amplitude/frequency limit (figure 2A, lower left). The most "stable" or commonly encountered phase-locking pattern was that of 1:1 labeled with a fuschia red color in figure 2A. Such low frequency firing patterns also have been reported for a broad range of current amplitudes when the driving stimulus was a sustained depolarization [21] (e.g., motoneurons). The Farey sequences for phase-locking patterns described above were conserved across the entire parameter space. The most prominent summation was the that of 2:1 (yellow) and 1:1 (fuschsia red) which yielded a region of 3:2 (persimmon orange). Smaller stable regions were observed interleaved between the larger regions. For example, the thin green region interleaved between the 3:2 and 1:1 regions corresponds to a 4:3 (mean nISI of 1.25; see color band). Another example of additivity is the thin blue 5:3 region between 3:2 and 2:1 phase-locking regions. Higher order phase-locking patterns, such as 11:1, 8:1, 6:1 and 4:1, were observed with the lower amplitude and frequency sinusoid driving. The latter regions were labeled in lavender on the parameter map of figure 2 (arrow). Higher order phase-locking patterns (e.g., see the lavender region) covered relatively less of the parameter space then their lower order counterparts (e.g., 1:1) as has been observed for analogous patterns generated by iterating simple (e.g., unimodal maps) mathematical functions [12,16]. Finally, with a low amplitude and high frequency of stimulus, multiple spikes occurred on each stimulus cycle. The latter patterns corresponded to fractional mean nISI's. Stable patterns of 1:2, 1:3, 1:4 and 1:5 with corresponding colors of pink, red, orange and yellow were observed.

Sinusoidally driven single neuron: aperiodic patterns

In the transition zones between stable phase-locking regions, there were regions with aperiodic entrainment patterns. The existence of such aperiodic regimes for the HH differential equations was predicted in a linearized version of these equations [3]; but the extent and nature of these aperiodicities had not been explored in detail. Standard deviations for aperiodic nISI sequences were higher than those associated with periodic sequences and they ranged from 1 to 6. Standard deviations of 0 were mapped in black and those with higher values were mapped on a rainbow color spectrum (figure 2B, color bar). Standard deviations greater than or equal to 6 were mapped in red. When larger steps in amplitude (1 ${mu A} over cm sup 2$) and $omega$ ($0.20 {msec sup -1}$) were used to map the frequency-amplitude parameter space, the highly variable region was located between $omega$ values of 0.70 to 0.85 ($msec sup -1$) or with spike frequencies of 111 to 135 Hz (figure 2B). However, if smaller steps in amplitude (0.5 ${nA} over cm sup 2$) were used, the aperiodic region extended over a much wider range of $omega$'s values (i.e., 0.21 and 0.80 $msec sup -1$) or frequencies (i.e., 33 135 Hz; figures 3 & 4). Thus, aperiodic patterns were observed within a broader range of stimulus frequencies then is immediately apparent from figure 2.

Aperiodic spiking patterns were very similar to the aperiodic patterns described for somatosensory fiber entrainment to sinusoid somatic stimulation [45]. In the first example, the aperiodic entrainment pattern was located between regions of 4:1 and 3:1 phase-locking for a stimulus period of 19.63 msec ($omega$= 0.33 $msec sup -1$; figure 3B, C). When the stimulus amplitude was increased by only 0.33 % (from 1504.0 to 1504.5 ${{nA} over cm sup 2}$) there was an associated increase in nISI mean (from 4 to 5.4) and standard deviation (from 0 to 1.5). The rise in standard deviation reflected the instability in spiking pattern. These patterns were considered "aperiodic" because no repeating sequence of ISI's was observed within the 10 sec steady-state period (figure 5, A & D). In the first example, spiking occurred at intervals ranging from 4 to 12 times the stimulus period. In the second example, the aperiodic pattern was located between stable phase-locking patterns of 11:1 and 8:1. In this case the stimulus period and amplitude were 20.94 ($omega$=0.34 $msec sup -1$) and 1540 (${{nA} over cm sup 2}$), respectively. Normalized ISI's in the second example ranged from 8 to 16 times the stimulus period. A 1 sec epoch of membrane voltage is plotted for each example of aperiodic entrainment to illustrate the membrane variability underlying the ISI variability (figures 3B and 4B).

For aperiodic entrainment patterns, the entire sequence of ISI's within the 10 sec steady-state epoch was not predictable; however, there were patterns in the sequential order of nISI's. The distribution of repeating or "near-repeating" short ISI sequences was examined by plotting each ISI against the ISI that followed in an ISI return map (ISI (i) versus ISI (i+1). Return maps revealed "near-repeating" nISI sequences (figure 5B, E). In the first example, there were extended sequences of nISI's near 4. Consequently, nISI's were heavily distributed about the point (4,4) on the return map (figure 5B). In contrast, there was never an incident where two or more nISI's of 5 occurred sequentially over the last 10 sec of a 50 sec simulation (figure 5B).

There were marked asymmetries in the densities and distributions of ISI's within an aperiodic sequence such that the most frequently observed ISI's were near those observed in stable entrainment patterns of neighboring regions of the parameter space. Asymmetries in the distribution densities of ISI's were readily visualized with density histogram plots (figure 5C, F). When stimulus amplitude and period were 1504.5 ($n A$) and 19.63 (msec) respectively, ISI's tended to be distributed about a limited set of near integer multiples of the stimulus period (e.g., 4, 5, 6, 12). In the first example, nISI's ranged from 4 to 12, with 4, 5 and 6 being the most frequent (largest peaks in the density histogram). A small decrease in stimulus amplitude was sufficient to generate a stable pattern of 4:1 or 5:1 (e.g., 4:1, figure 3A). Thus, the most frequently observed ISI's in aperiodic patterns tended to be of similar normalized duration to those observed in the surrounding regions of parameter space. In the second example, nISI's were distributed about 9, 10, and 11 most frequently (figure 5F, right). Similarly, a stable pattern of 11:1 was observed in neighboring regions of the parameter space (figure 4A). In summary, even though ISI sequences were unpredictable for a given aperiodic entrainment pattern, the range and density of ISI's were constrained by the stimulus parameters used.

Entrainment Patterns for a Collection of Neurons

What happens to the dynamics of a sinusoidally driven neuron if it is embedded within a matrix of synaptically coupled neurons? A network of 25 coupled neurons was constructed to explore this problem. Neurons were modeled with the HH equations as above. Each neuron received synaptic input from 5 to 10 of its nearest neighboring neurons via a rapid depolarizing postsynaptic current (PSC, see Methods). The model PSC followed each action potential with a short delay which was 3 to 15 msec, depending upon displacement between input and target neuron. A 3 msec delay between neighboring model neurons approximated experimentally observed intracortical and interhemispheric monosynaptic delays [46] and allowed for a rich array of autonomous firing patterns (Siegel and Read, unpublished data). These neurons lacked the multiple compartments and sites of PSC input which can create additional temporal delays. The given parameters were selected so that this small highly interconnected network was capable of generating a sustained output after one or more of the 25 neurons was stimulated to threshold. The range and degree of synchronized or asynchronized activity of the population model did not change appreciably when the number of neurons was increased to 100; hence, the 25 neuron model was used. Once threshold was reached for the central neuron, there was a short period of time (up to 1 sec) before a spiking pattern was considered stable (e.g., figure 6, bottom). Thus, only those ISI's occurring in the last 9 sec of a 10 sec simulation were used to assess entrainment periodicity in the population model.

Neurons embedded in the network fired either phase-locked to the sinusoid rhythm, or in an aperiodic fashion, depending upon the sinusoid period and amplitude. The amplitude of the sinusoid stimulus to cell 13 was varied from 1550 to 1600 (${{nA} over cm sup 2}$) in 0.5 nA increments. Entrainment to the sinusoid rhythm was dependent upon the stimulus amplitude and frequency (e.g., figure 6, top). When the embedded neuron was periodically entrained to the driving stimulus, nISI's were equal to the nISI mean and the standard deviation was zero. When the neuron failed to periodically entrain, nISI's were scattered over a range of values and were associated with an increased standard deviation (figure 6, open circles). Sinusoid amplitudes and frequencies which resulted in phase-locked patterns for cell 13 also entrained the entire population to a common rhythm. Similarly, those stimulus conditions that generated aperiodic output in the driven neuron resulted in an asynchronous mode of output from the population. For example, the aperiodic activity of neuron 13 in figure 6 (bottom) was associated with asynchronous output from the coupled population (not shown). Hence, the entrainment pattern of a single neuron (cell 13) reflected the entrainment activity of the entire population.

Embedding the sinusoid driven neuron within an excitatory network dramatically restricted the range of entrainment patterns. Higher order phase-locked patterns (e.g., low frequency entrainment such as 11:1, 5:1, 3:1) were no longer observed in the population model. The most prominent phase-locked pattern in the population coupled neuron was a pattern of 1:1 phase locking with nISI's near a mean of 1. When neuron 13 was driven with a sinusoid period of 19.6 msec, aperiodic ISI sequences were common. The ISI's within an aperiodic sequence were restricted to a range of 3.92 to 39.2 msec or a nISI range of 0.2 to 2 (sinusoid period = 19.6 msec; figure 6, top). It is noteworthy that this is a biologically feasible range of entrainment frequencies, as ISI's of 39.2 to 3.92 msec correspond to frequencies of 25 to 255 Hz. The short ISI's associated with periodic and aperiodic entrainment in the population were a reflection of the short delay in feedback from neighboring neurons.

Discussion

Sinusoidal stimulated HH differential equations were used to model the basic features of entrainment of an isolated neuron and a neuron embedded within a population of synaptically coupled neurons. The resting state of the HH neuron was quiescent, that is, there was no baseline spiking prior to sinusoid current input. ISI's were normalized to the stimulus period, and means and standard deviations of nISI's were used as measures of entrainment. Spiking of HH neurons was phase-locked for most stimulus conditions. Systematic variation of the stimulus amplitude and period revealed ordered transitions between mean ISI's and phase-locked patterns, the order of which was predictable according to a simple additive rule (Farey summation). Higher order and aperiodic entrainment patterns were observed interleaved between the more stable phase-locking domains of the parameter space. Thus the hypothesis that aperiodic entrainment patterns arise in part from the membrane conductance dynamics was supported. However, aperiodic patterns of entrainment were more prevalent when there was feedback from a population affecting the periodically driven neuron.

Farey summation between neighboring phase-locking domains within the frequency-amplitude parameter space is a general property of both "threshold" and "coupled" oscillatory systems [1,3,16,27]. $Threshold$ oscillatory systems are those systems that have a threshold for generating oscillatory output [3,27]. Many nerve fibers or neurons are quiescent in the resting state and can be considered to have a threshold for oscillatory (spike) output [4,11,33,40,45]. Systems that oscillate intrinsically at rest are considered $coupled$ oscillatory systems when driven with an extrinsic rhythmic input. A resting or baseline oscillatory state is the natural state of cardiac fibers [16] and it is experimentally induced in neurons by lowering the external calcium concentrations [1,27] or by injecting depolarizing bias current [8,19,21]. In both systems the phase-locked entrainment patterns vary systematically with changes in stimulus amplitude and frequency [1,17,19,21,27,31,32] and these entrainment patterns can be predicted with the simple Farey summation rule [12,15,39] (Results and figure 2).

The frequency-amplitude parameter space of the $threshold$ system is sufficiently different from that of $coupled$ oscillatory systems so that the parameter space effectively serves as a "fingerprint" of the system. In $threshold$ systems, such as the quiescent HH neuron, input currents of low frequency and low amplitude fail to evoke sodium spikes or resulting ISI's (figure 2A). As either frequency or amplitude is systematically increased there is a monotonic increase in firing rate. Neighboring regions of different entrainment patterns originate and diverge from a common (or critical) point in the frequency-amplitude parameter space. Neighboring regions with different entrainment patterns abut one another so that changes in entrainment pattern with stimulus parameters are abrupt. Conversely, when an HH neuron (or cardiac fiber) is placed into a state of rhythmic activity prior to entrainment, the boundaries of entrainment regions are less discrete. A systematic increase in frequency and amplitude does not lead to a monotonic increase in firing rate. In $coupled$ oscillatory systems, the different entrainment regions no longer arise from a common point; they instead arise from separate points to form "tongue" like structures [1,15,17,32]. Transitions from one entrainment region to the next are less abrupt. This separation of entrainment regimes is a hallmark difference between the two types of systems. We refer to the convergence of the $threshold$ parameter space into tightly packed phase-locking regions as a "folding" of the parameter space and conversely to the separation of these regimes in the $coupled$ oscillatory system as an "unfolding".

Aperiodic entrainment regimes have been well documented for entrainment of $coupled$ oscillatory systems [17,19,21,32,38,39]; however, they are less well documented for the $threshold$ systems. Aperiodic entrainment regimes were observed in the folded parameter space of the HH neuron with a threshold for output (figures 2-4); thereby, confirming predictions from 1) experimental systems and 2) a simple mathematical model. Systematic transitions from periodic to aperiodic firing have been demonstrated for neurons that fire infrequently or not at all prior to stimulus input [19,21,31,45]. Furthermore, the genesis of aperiodic patterns from the HH differential equations (with a threshold) had been predicted from a linearized model, the Fitzhugh-Nagumo equations [3]. That the higher order phase-locking and interleaved aperiodic regimes were conserved in the more realistic HH model suggests that the origin of aperiodicities from the system is a robust phenomenon (i.e., these regions were not obscured by small integration errors). The exact temporal sequence for aperiodic output, given a set of parameters, was not predictable. However, the frequency distribution of ISI's for aperiodic regimes overlapped with those of neighboring phase-locked regions of frequency-amplitude parameter space (e.g., figures 3 and 4). Hence, the mean frequency of spiking in aperiodic regimes was somewhat predictable. Larger aperiodic regimes were observed on the borders between phase-locking domains with HH equations (and simpler models) that have been adjusted to allow for sustained oscillation prior to periodic stimulation [17,19,21,32,38,39]. It had been suggested that the presence of a "carrier" frequency or baseline oscillatory state increased these regions of aperiodicity, though a mechanism was not offered [21]. Since aperiodic patterns in either system tended to occur on the borders of phase-locking domains, it is possible that the increase in aperiodic domains stemmed from the unfolding of the frequency-amplitude parameter space observed in $coupled$ oscillatory systems.

If the presence of an intrinsic "carrier" frequency or baseline output of a single neuron has the effect of unmasking new aperiodic entrainment patterns, then what effect does feedback from a population of neurons have? Unlike the primary somatosensory fibers, most secondary sensory or central neurons are synaptically $coupled$ to a matrix of neurons. In the present model, a single neuron was driven by sinusoid input to emulate a periodic afferent drive (or sensory stimulus). Embedding the sinusoidally driven neuron within a population of neurons dramatically altered the model neuron's entrainment properties and increased the incidence of aperiodic entrainment. High order, low frequency phase-locking was eliminated and regions of aperiodicity were expanded in the frequency-amplitude parameter space (compare figures 1 and 6). Furthermore, the output during aperiodic or periodic entrainment was restricted to a frequency range of 25 to 255 Hz. Thus, in this simplified population model aperiodic entrainment was a dominant pattern and spiking fell within a narrow band of frequencies.

Aperiodic or 'noisy' entrainment to a periodic stimulus is a common feature of $in$ $vivo$ neural ensembles; however the source and nature of such aperiodicities has not been determined. Rhythmic synaptic drive or membrane potential oscillations are ubiquitous phenomena in nervous systems [2,6,9,14,24,45]; and yet neuronal spike train output is generally not periodic itself, i.e., entrainment is aperiodic. In the somatosensory system, aperiodic entrainment is observed in primary sensory fibers and is even more prevalent in second order sensory neurons [29,45]. A subpopulation of hippocampal and entorhinal neurons fires in phase with a 200 Hz afferent drive (ripple) [6]. The output of hippocampal neurons is correlated but aperiodic or intermittent with respect to the 200 Hz afferent drive. Clearly one source of aperiodicity in such neuronal populations is the random membrane fluctuations caused by synaptic or membrane channel "noise" [7,26,42](c.f., Liebovitch, 1993). Accordingly, spontaneous interspike interval or synaptic potential data from some sensory and motor neurons can be approximated by Poisson or Gaussian distributions [13,35,37]; and the data distribution can be approximated in simulations by adding a stochastic element to account for random membrane fluctuations [7,19,26,41,42,44]. However, simply driving a model or real neuron with a random input is not a sufficient way to generate highly variable spike trains under all stimulus conditions (e.g., low amplitude random current fluctuations [35,44]). Indeed, evidence to the contrary exists [21,44,47] and suggests that either 1) Poisson fits of experimental spike data have to be questioned [47] or 2) existing models for simulating such distributions have to be modified. A second source of $jitter$ for entrainment of sensory neurons could be non-stochastic properties such as those modeled here including: 1) inherent membrane properties, 2) synaptic potential kinetics and 3) axonal conductance delays [13,17,45]. It will be of great interest in future studies to 1) identify unique frequency ranges of biological $jitter$ and 2) determine the role, if any, such $jitter$ has for information processing in the brain.

In the past, biologists have studied physiologic mechanisms for improving signal to noise ratios with the tacit assumption that we, as observers, know the difference between signal and noise in the brain [35]. The change in mean firing rate that occurs immediately following sensory stimulation is generally considered a meaningful signal in sensory neurons at all levels including cortical. Hence, manipulations which decrease "baseline" firing rates and increase stimulus specific firing rates are categorized as signal enhancing [28]. Unfortunately, stimulus specific changes in mean firing are operationally defined and other 1) sensory event and 2) neural network related spiking patterns may go undetected when mean firing rate is the sole measure of signal [35]. Accordingly, Vaadia and colleagues [48] found behavior related changes in temporal patterns which were not associated with changes in mean firing rate and Nakamura [30], Optican [34] and their respective colleagues found stimulus specific temporal spike or burst patterns. Thus, distinct stimuli for sensory neuron activation may be associated with emergence of unique temporal patterns of output. Spike train patterns need not be transmitted spike for spike in order for temporal patterns to be meaningful to target neurons. Indeed, higher order spiking patterns generated when a neuron (or population of neurons) fails to spike on every stimulus cycle (e.g., 3:2, or aperiodic) may more effectively activate, or conceivably potentiate, recipient neurons [22,23,35].

In conclusion we have shown that there are differences in the patterns of activity that may evolve from individual neurons and populations of interconnected neurons. The single neuron tends to follow the stimulus in a periodic phase-locked pattern, whereas, the population tends to fire aperiodically over a large stimulus regime. Such irregular dynamics are often found $in$ $vivo$ in neurons that are part of highly interconnected neural networks. However, it remains an open question whether these complicated dynamics recorded in the brain arise in part from deterministic interactive mechanisms, as posited here, or solely from the presence of stochastic processes.


References

1. Aihara, K. (1984) Periodic and non-periodic responses of a periodically forced Hodgkin-Huxley oscillator. J. Theor. Biol. 109: 249-269.

2. Alonso, A., and Klink, R. (1993) Differential electroresponsiveness of stellate and pyramidal-like cells of medial entorhinal cortex layer II. J. Neurophys. 70: 128-143.

3. Alexander, J.C., Doedel, E.J., and Othmer, H.G. (1990) On the resonance structure in a forced excitable system. SIAM J. Appl.Math. 50: 1373-1418.

4. Bishop, P.O. (1984) Processing of visual information within the retinostriate system. In: The handbook of physiology: The nervous system III (section 1). Waverly Press, Inc., Baltimore.

5. Brown, T.H., and Johnston, D. (1983) Voltage-clamp analysis of mossy fiber synaptic input to hippocampal neurons. J. Neurophysiol. 50: 487-507.

6. Buzsaki, G., Bragin, A., Chrobak, J.J., Nadasdy, Z., Sik, A., Hsu, M., and Ylinen, A. (1994) Oscillatory and intermittent synchrony in the hippocampus: relevance to memory trace formation. In: Temporal Coding in the Brain, Springer-Verlag Berlin.

7. Calvin, W.H., and Stevens, C.F. (1968) Synaptic noise and other sources of randomness in motoneuron interspike intervals. J. Neurophysiol. 31, 547-587.

8. Canavier, C.C., Clark, J.W., and Byrne, J.H. (1990) Routes to chaos in a model of a bursting neuron. Biophys. J. 57: 1245-1251.

9. Carr, C.E., Heiligenberg, W., and Rose, G.J. (1986) A time-comparison circuit in the electric fish midbrain. I. behavior and physiology. J. Comp. Neurol. 314:306-318.

10. Coon, D.D., Ma, S.N., and Perera, G.U. (1987) Farey-fraction frequency modulation in the neuronlike output of silicon p-i-n Diodes at 4.2 K. Physical Rev. Lett. 58: 1139-1142.

11. Darian-Smith, I. (1984) Thermal sensibility. In: The handbook of physiology: The nervous system III (section 2). Waverly Press, Inc., Baltimore.

12. Feigenbaum, M. (1988) Complicated objects on regular trees. In: Nonlinear evolution and chaotic phenomena. Plenum Publishing Corp.

13. Freeman, W.J. (1975) Mass action in the nervous system. New York: Academic Press.

14. Fregnac, Y., Bringuier, V. and Baranyi, A. (1994) Oscillatory neuronal activity in visual cortex: a critical re-evaluation. In: Temporal Coding in the Brain, Springer-Verlag Berlin.

15. Glass, L., and Mackey, M.C. (1988) Periodic stimulation of biological oscillators. In: From clocks to chaos. The rhythms of life., Princeton University Press.

16. Glass, L., Guevara, M.R., and Shrier, A. (1987) Universal bifurcations and the classification of cardiac arrhythmias. Ann. N.Y. Acad. Sci. 504: 168-178.

17. Guevara, M.R., Glass, L., Mackey, M.C., and Shier, A. (1983) Chaos in neurobiology. IEEE Trans on Sys, Man, Cybern 5: 790-798.

18. Guevara, M.R., Glass, and Shier, A. (1981) Phase-locking, period-doubling bifurcations, and irregular dynamics in periodically stimulated cardiac cells. Science 214: 1350-1353.

19. Guttman, R., Feldman, L., and Jakobsson, E. (1980) Frequency entrainment of squid axon membrane. J. Membrane Biol. 56: 9-18.

20. Hodgkin, A.L., and Huxley, A.F. (1952) A quantitative description of membrane current and its application to conduction and excitation in nerve. J. Physiol. 117: 500-544.

21. Holden, A.V. (1976) The response of excitable membrane models to a cyclic input. J. Cybernetics 21: 1-7.

21a. Jack, J.J.B., Noble, D., and Tsien, R.W. (1975) Electric Current Flow in Excitable Cells. London: Oxford Univ. Press.

22. Klemm, W.R., and Sherry, C.J. (1982) Do neurons process information by relative intervals in spike trains? Neurosci. & Biobehav. Rev. 6: 429-437.

23. Larson, J. and Lynch, G. (1986) Induction of synaptic potentiation in hippocampus by patterned stimulation involves two events. Science 232: 985-988.

24. Llinas, R.R., Grace, A.A., and Yarom, Y. (1991) In vitro neurons in mammalian cortical layer 4 exhibit intrinsic oscillatory activity in the 10- to 50-Hz frequency range. Proc. Natl. Acad. Sci. 88: 897-901.

25. Liebovitch, L.S. (1993) Interpretation of protein structure and dynamics from the statistics of the open and closed times measured in a single ion-channel protein. J. Stat. Physics 70: 329-337.

26. Longtin, A. (1993) Stochastic resonance in neuron models. J. of Stat. Phys. 70:309-327.

27. Matsumoto, G., Aihara, K., Hanyu, Y., Takahashi, N., Yoshizawa, S. and Nagumo, J. (1987) Chaos and phase locking in normal squid axons. Physics A 123: 162-165.

28. McCormick, D.A. (1989) Cholinergic and noradrenergic modulation of thalamocortical processing. TINS 12: 215-221.

29. Mountcastle, V.B., Steinmetz, M.A., and Romo, R. (1990) Frequency discrimination in the sense of flutter: psychophysical measurements correlated with postcentral events in behaving monkeys. J. Neurosci. 10: 3032-3044.

30. Nakamura, K., Mikami, A. and Kubota, K. (1991) Unique oscillatory activity related to visual processing in the temporal pole of monkeys. Neurosci. Res. 12: 293-299.

31. Nemoto, I., Miyazaki, S., Saito, M. and Utsunomiya, T. (1975) Behavior of solutions of the Hodgkin-Huxley equations and its relation to properties of mechanoreceptors. Biophysical J. 15: 469-479.

32. Nomura, T., Sato, S., Doi, S., Segundo, J.P. and Stiber, M.D. (1993) A Bonhoeffer-van Pol oscillator model of locked and non-locked behaviors of living pacemaker neurons. Biological Cybernetics 69: 429-437.

33. Norgren, R. (1984) Central neural mechanisms of taste. In: The handbook of physiology: The nervous system III (section 2). Waverly Press, Inc., Baltimore.

34. Optican, L.M., and Richmond, B.J. (1987) Temporal encoding of two-dimensional patterns by single units in primate inferior temporal cortex, J. Neurophysiol. 57: 167-178.

35. Perkel, D.H. and Bullock, T.H. (1968) Neuronal Coding: a report based on an NRP work session. Neurosciences Res. Prog. Bull. 6: 221-348.

36. Press, W.H., Flannery, B.P., Teukolsky, S.A. and Vetterling, W.T. (1988) Integration of ordinary differential equations. In: Numerical recipes. Cambridge University Press, Cambridge.

37. Prucnal, P.R. and Teich, M.C., (1983) Refractory effects in neural counting processes with exponentially decaying rates. IEEE Trans on Sys, Man, and Cybern

38. Rinzel, J. and Miller, R.N. (1980) Numerical calculation of stable and unstable periodic solutions to the Hodgkin-Huxley equations. Math. Biosci., 49: 27-59.

39. Sato, S. and Doi, S. (1992) Response characteristics of the BVP neuron model to periodic pulse inputs. Mathematical Biosciences, 112: 243-259.

40. Scheibel (1984) The brain stem reticular core and sensory function. In: The handbook of physiology: The nervous system III (section 1). Waverly Press, Inc., Baltimore.

41. Siegel, R.M. (1990) Non-linear dynamical system theory and primary visual cortical processing. Physica D 42: 385-395.

42. Siegel, R.M. and Read, H.L. (1993) Models of the temporal dynamics of visual processing. J. Stat. Physics 70: 297-308.

43. Siegel, R.M., Tresser, C. and Zettler, G. (1992) A decoding problem in dynamics and in number theory. Chaos 2: 473-494.

44. Stein, R.B. (1967) Some models of neuronal variability. Biophysical Journal 7: 37-68

45. Talbot, W.H., Darian-Smith, I., Kornhuber, HH, and Mountcastle, VB (1968) The sense of flutter-vibration: comparison of human capacity with response patterns of mechanoreceptive afferents from the monkey's hand. J. Neurophysiol. 31: 302-334.

46. Thomson, A.M. (1986) A magnesium-sensitive post-synaptic potential in rat cerebral cortex resembles neuronal responses to N-methylaspartate. J. Physiol. 370: 531-549.

47. Teich, M.C. and Khanna, S.M. (1985) Pulse-number distribution for the neuronal spike train in the cat's auditory nerve. J. Acoust. Soc. Am. 77: 1110-1128.

48. Vaadia, E., Haalman, I., Abeles, M., Bergman, H., Prut, Y., Slovin, H., and Aertsen, A. (1995) Dynamics of neuronal interactions in monkey cortex in relation to behavioral events. Nature 373: 515-518.