The Journal of Neurophysiology Vol. 77 No. 5 May 1997, pp. 2736-2752
Copyright ©1997 by the American Physiological Society

Low-Amplitude Oscillations in the Inferior Olive: A Model Based on Electrical Coupling of Neurons With Heterogeneous Channel Densities

Yair Manor^{1, 2}, John Rinzel³, Idan Segev^{1, 2}, and Yosef Yarom^{1, 2}

¹ Neurobiology Department, Hebrew University, Jerusalem 91904, Israel; ² Interdisciplinary Center for Neural Computation, Hebrew University, Jerusalem 91904, Israel; and ³ National Institute of Diabetes and Digestive and Kidney Diseases, Bethesda, Maryland 20814

	ABSTRACT

Abstract Introduction Methods Results Discussion References

Manor, Yair, John Rinzel, Idan Segev, and Yosef Yarom. Low-amplitude oscillations in the inferior olive: a model based on electrical coupling of neurons with heterogeneous channel densities. J. Neurophysiol. 77: 2736-2752, 1997. The mechanism underlying subthreshold oscillations in inferior olivary cells is not known. To study this question, we developed a single-compartment, two-variable, Hodgkin-Huxley-like model for inferior olive neurons. The model consists of a leakage current and a low-threshold calcium current, whose kinetics were experimentally measured in slices. Depending on the maximal calcium and leak conductances, we found that a neuron model's response to current injection could be of four qualitatively different types: always stable, spontaneously oscillating, oscillating with injection of current, and bistable with injection of current. By the use of phase plane techniques, numerical integration, and bifurcation analysis, we subdivided the two-parameter space of channel densities into four regions corresponding to these behavioral types. We further developed, with the use of such techniques, an empirical rule of thumb that characterizes whether two cells when coupled electrically can generate sustained, synchronized oscillations like those observed in inferior olivary cells in slices, of low amplitude (0.1-10 mV) in the frequency range 4-10 Hz. We found that it is not necessary for either cell to be a spontaneous oscillator to obtain a sustained oscillation. On the other hand, two spontaneous oscillators always form an oscillating network when electrically coupled with any arbitrary coupling conductance. In the case of an oscillating pair of electrically coupled nonidentical cells, the coupling current varies periodically and is nonzero even for very large coupling values. The coupling current acts as an equalizing current to reconcile the differences between the two cells' ionic currents. It transiently depolarizes one cell and/or hyperpolarizes the other cell to obtain the regenerative response(s) required for the synchronized oscillation. We suggest that the subthreshold oscillations observed in the inferior olive can emerge from the electrical coupling between neurons with different channel densities, even if the inferior olive nucleus contains no or just a small proportion of spontaneously oscillating neurons.

	INTRODUCTION

Abstract Introduction Methods Results Discussion References

The olivocerebellar system, which is involved in motor control (Holmes 1939; Llinás 1984; Llinás and Welsh 1993) and may also participate in motor learning (Albus 1971; Marr 1969; Robinson 1976), generates a rhythmic activity at a frequency of ~10 Hz. This rhythmic activity is expressed as a temporal relationship of complex spike activity in the cerebellum (Llinás and Sasaki 1989), and, under pathological conditions, it is manifested as an "enhanced physiological tremor" (Llinás 1984). It has been postulated that the subthreshold spontaneous (sinusoidal-like) oscillations (STOs) of membrane potential in inferior olivary (IO) neurons, observed in slice preparations, underlie this rhythm (Bloedel and Ebner 1984; Llinás and Sasaki 1989). The STO in an IO neuron acts to rhythmically change its firing probability. Consequently, the target cerebellar Purkinje cells are activated, with some probability, in particular time windows. Not every Purkinje cell will definitely receive adequate input to fire on each cycle, however. Thus the rhythmic activity and its relationship to motor behavior can be revealed only after the activity is recorded simultaneously from a large number of Purkinje neurons (Llinás and Sasaki 1989; Welsh et al. 1995). Although attempts to characterize this rhythmicity have been made, they were unsuccessful when single units were recorded (Keating and Thach 1995). In addition to their proposed role in generating the cerebellar rhythmic activity, the STOs may also act to synchronize the inputs to the cerebellum under normal conditions (Lampl and Yarom 1993).

Several observations are noteworthy. The STO frequency ranges over 4-10 Hz, and the amplitude varies between 3 and 10 mV (Llinás and Yarom 1986a). STOs were observed in only 10% of the slice preparations; in each oscillating slice, they were recorded in most of the neurons encountered (Lampl and Yarom 1997; Yarom 1991). They are sensitive to calcium blockers, but are unaffected by application of sodium blockers (Benardo and Foster 1986; Llinás and Yarom 1986a). Octanol, a low-threshold calcium current blocker (Llinás and Yarom 1986b), blocks the STOs (Lampl and Yarom 1997). Gross extracellular stimulation affects the oscillations, but intracellular stimulation of any given neuron does not (Llinás and Yarom 1986a). On the other hand, global hyperpolarization of the neurons (by decreasing the extracellular K⁺ concentration) reversibly blocks the STO (Lampl and Yarom 1996). In cases where the membrane potential is not oscillating, neither extracellular nor intracellular stimuli can make the impaled neuron oscillate. Yet, intracellular injection of a hyperpolarizing current pulse may generate a rebound response on release to rest. This rebound response consists of one or more low-threshold calcium spikes. The amplitudes of these responses decrease successively; the neuron behaves like a damped oscillator (Yarom 1991).

Because IO cells are electrically coupled (de Zeeuw et al. 1990; Llinás and Yarom 1981; Llinás et al. 1974; Sotelo et al. 1974), it has been hypothesized that the IO network is composed of damped oscillators that, when coupled, generate a sustained oscillatory behavior (Yarom 1991). Several lines of evidence support this hypothesis. First, developmental studies show a clear correlation between the times of gap junction formation and the development of STOs in rats (Bleasel and Pettigrew 1992). Second, simultaneous recordings from two neurons show that nearby neurons oscillate with the same frequencies and phases (Benardo and Foster 1986; Llinás and Yarom 1986a). Third, pharmacological treatments that block electrotonic transmission, such as exposure to bicarbonate-buffered solution, abolish the STO (Bleasel and Pettigrew 1994). Fourth, even if some individual cells oscillate spontaneously, there is no evidence that the proportion is large enough that these cells could act as a pacemaker to drive a network rhythm.

Isopotential quiescent cells that are identical in all properties cannot mediate sustained synchronized oscillations when coupled via linear-resistive gap junctions. Thus we have developed a heterogeneity hypothesis that could explain the generation of STOs in the IO. We constructed an experimentally based yet minimal model of IO cells that contains only two currents: a low-threshold, inactivating calcium current and a passive leakage current. We showed that different combinations (amounts) of these two currents can support markedly distinct electroresponsiveness of the IO cells, including cells with a unique stable resting state (possibly with damped oscillatory transients), spontaneously oscillating cells, or bistable cells with plateauing behavior. Utilizing the model, we demonstrate how two different nonoscillating neurons can generate low-amplitude oscillations when they are electrically coupled.

	METHODS

Abstract Introduction Methods Results Discussion References

Extraction of voltage-clamp data

The low-threshold calcium current typical of IO neurons is the major conductance responsible for the generation of STOs (Lampl and Yarom 1996; Llinás and Yarom 1986a; Manor 1995). We extracted the gating kinetics of this conductance by voltage-clamp experiments. A detailed description of the experimental results, the space-clamp problems, experimental protocols, and theoretical validations is given elsewhere (Manor 1995). With voltage-clamp data from 15 cells, we determined an activation curve (m³) and an inactivation curve (h) as shown in Fig. 1A. The solid curves are the equations

(1)

and

(2)

The similarities between the results for the different cells further support the accuracy of our voltage-clamp protocols. The dependence of the inactivation time constant (_h) on voltage was extracted from six cells (Fig. 1B). It is comparable with the time constant of inactivation found in thalamic relay neurons (Coulter et al. 1989; Huguenard and Prince 1992). Our data are best fit (solid curve) with the following bell-shaped equation

(3)

View larger version (17K): [in this window] [in a new window]   FIG. 1. Experimental extraction of activation and inactivation curves. A: activation (m3, , , , , ) and inactivation (h, , , , , ) data from several different olivary neurons. Solid curves: Boltzman functions (Eq. 1 and 2) that are fitted to these data. Experimental details are given in Manor (1995). B: data for the time constants of inactivation of several olivary neurons. Solid curve: best fit to the data (Eq. 3).

The activation time constant, _m, was found to range between 5 and 15 ms, an order of magnitude faster than _h. Thus we approximated the activation as an instantaneous function of the membrane potential. That is, m  m(V).

SINGLE-NEURON MODEL. Two differential equations govern the dynamics of the single cell

(4)

and

(5)

where V is the membrane potential (in mV), C_m = 1 µF/cm² is the specific capacitance of the membrane, and  is the temperature correction factor (which we set to 1), I_app is the applied current, and I_ion is the sum of two ionic currents, a low-threshold calcium current I_T and a leakage current I_L (all currents in µA/cm²)

(6)

I_T is described by the Hodgkin-Huxley formalism with two gating variables: the rapid activation variable m and the inactivation variable h

(7)

where _T is the maximal calcium conductance (in mS/cm²) and V_Ca is the calcium reversal potential, set to +120 mV. The leakage current I_L is modeled by where g_L is the leak conductance and V_L is the reversal potential of the leakage current, set to 63 mV.

Tools for analyzing the model's dynamic behavior

We are interested in studying the oscillatory (periodic) and steady-state (time-independent) modes of behavior of our neuron model as functions of the channel densities and injected current. As parameters are varied, the behavior may change from one mode to another if a solution loses stability; a system subject to physiological random noise cannot reside for long in an unstable state. In some regimes, multiple modes coexist and, depending on the initial conditions (the initial values of V and h), the system may converge to any one of the stable (i.e., attracting) modes. Thus it is important to determine the stability of the solution. Here we briefly indicate aspects of considering the stability of steady-state solutions. By definition of a steady state, each variable's time derivative equals zero. One asks whether a small perturbation from this steady state will decay to zero with time (the system returns to the steady state) or grow, leading the system into a different behavior. In the former case, the steady state is stable; in the latter it is unstable.

According to linearized perturbation theory, for a two-variable cell model the stability is determined by the eigenvalues of the two-dimensional Jacobian matrix. If the eigenvalues have negative real part, the steady state is stable; it is unstable if either eigenvalue has positive real part. When the eigenvalues have an imaginary part, time-dependent trajectories near the steady state are oscillatory, either damped or growing. As parameters are modified, if the real part of a complex pair of eigenvalues changes from negative to positive, the steady state changes from stable to unstable. Moreover, at the transition point the eigenvalues sum to zero. Thus, because in general the eigenvalue sum equals the sum of the Jacobian's diagonal elements, the following condition is satisfied by our system at such a transition point

(8)

where I_ion/V is the slope of the instantaneous current-voltage (I-V) relation (Rinzel and Ermentrout 1989, their Eq. 5.19). This point is called a Hopf bifurcation (HB). In general, a periodic solution, also called a limit cycle, emerges. The oscillation's frequency equals 2 times the imaginary part of the eigenvalues; the emergent oscillation may be stable or unstable, depending on parameter values. For a more detailed description of linear perturbation theory methods, see Edelstein-Keshet (1988) and Strogatz (1994).

The steady states and their stability for our system were also explored graphically by phase plane methods (see Edelstein-Keshet 1988 and Strogatz 1994 for general background). A phase plane portrait is a two-dimensional plot describing the dynamic relationships between the two dependent variables, in our case V and h;solution trajectories are curves in the V-h plane. The nullclines of the system are the two special curves (not trajectories) along which

(9)

From Eq. 5, the h nullcline is the curve

(10)

Denoting the relationship between V and h along the V nullcline by h = _V(V), we see that it is defined implicitly from Eq. 4

(11)

Important insights may be obtained by plotting the nullclines and examining their relationship (Fitzhugh 1969; Rinzel and Ermentrout 1989). Phase planes in this paper are found in Fig. 5; here, one sees the usual N-shaped V nullcline of an excitable system. An intersection of the two nullclines defines a steady-state point (, ). At critical parameter values where an HB occurs, Eq. 8 is satisfied. Thus a necessary condition for Hopf-like instability of (, ) is

(12)

(Rinzel and Ermentrout 1989). To interpret this condition geometrically, we rewrite it in terms of the V nullcline's slope. Differentiating Eq. 11 yields an expression for the slope _V where  denotes derivative with respect to V. Thus, with the use of the chain rule for differentiation, we get This equation relates the slope of the I-V relation to the V nullcline's slope

(13)

In the present model In particular, applying these expressions at the steady state (, ), we can express the instability condition (Eq. 12) in the following equivalent way

(14)

where _V is the slope of the V nullcline (h = _V) at .

View larger version (14K): [in this window] [in a new window]   FIG. 5. Similarity between infinite and strong coupling. The dynamics of two strongly coupled cells are compared with the behavior of the average cell. Cells 1 and 2 are as in Fig. 4A (T = 0.4 mS/cm2, gL = 0.1 mS/cm2 and T = 0.4 mS/cm2, gL = 0.2 mS/cm2, respectively). A: V nullclines of cell 1 (· · ·), cell 2 (- - -), and the average cell (T = 0.4 mS/cm2, gL = 0.15 mS/cm2; ). The average cell is electrically equivalent to the case of cell 1 coupled to cell 2 with infinite coupling. In all 3 cases, the h nullcline () is identical. In each case, the equilibrium (fixed point) occurs at the intersection of the V nullcline with the h nullcline (59.8, 52.8, and 56.6 mV for cells 1 and 2 and average cell, respectively). In the cases of cell 1 and cell 2, the slope of the V nullcline at the intersection is positive. In the case of the average cell, the slope is negative and, with the parameters of the model, it satisfies the condition for instability (Eq. 14). B and C: cell 1 and cell 2 are electrically coupled with gcoup = 0.5 mS/cm2. In B, the limit cycle trajectories of cell 1 (· · ·) and cell 2 (- - -) and the limit cycle trajectory of the average cell () are superimposed on the nullclines of the average cell (). Arrows on limit cycle trajectories: direction of motion. In C, the ionic currents of cell 1 (· · ·) and cell 2 (- - -) and the coupling currents (: Icoup, from cell 2 to cell 1, and Icoup, from cell 1 to cell 2) are plotted as functions of time.

In other words, a steady state (for such a two-variable model) can be unstable only if it occurs where the V nullcline has a slope that is sufficiently negative; or equivalently, from Eq. 12, where the membrane's instantaneous I-V relation has sizable negative resistance. Moreover, destabilization cannot occur if the rate of inactivation (C_m/_h) is too large or the maximal calcium conductance (_T) is too small (see Eq. 14). By considering how the nullclines change with parameter values, one can identify ranges where multiple steady states or oscillations may or may not exist. Also, examination of the phase plane and nullclines allows one to know where V and h are increasing or decreasing and thereby to approximately predict the form of trajectories, such as limit cycles.

We used the software XPPAUT by B. Ermentrout to simulate the time courses of V and h and to carry out phase plane analysis by computing nullclines and stability of steady states. XPPAUT (available at ftp://ftp.math.pitt.edu/pub/bardware) also incorporates important features of the AUTO package (Doedel 1981). AUTO can perform bifurcation analysis of systems of ordinary differential equations. Among its capabilities are tracing stationary (steady-state) solutions as a parameter of the model is continuously changed, detecting limit points (such as HB points, turning points, period doubling bifurcations, bifurcations to tori, etc.), tracing periodic solutions, and following a limit point as two parameters of the model are free to change (a two-parameter continuation problem). We used AUTO to explore the behavior of an isolated cell (a two-variable model), or a pair of coupled cells (a four-variable model), when different parameters of the models are modified. In the latter case, the stability of the full system is estimated by computing the eigenvalues of a 4 × 4 Jacobian matrix. We focused on changes in _T, g_L, and I_app because these were the only parameters we modified from run to run, and from cell to cell in the same run.

Transient behavior in voltage-clamp mode

For an additional viewpoint on the dynamic stability of membrane potential steady states, more familiar to some physiologists, we used NEURON (Hines 1993) to simulate voltage- and current-clamp experiments and to construct momentary voltage-current (V-I) relationships, as defined in Jack et al. (1983). Resting values of the membrane potential and the inactivation variable were determined by an iterative procedure. Then the voltage was stepped from the resting state to some clamped value. At some specified time t after the onset of the voltage step, we measured the total ionic current (I_T + I_L). The momentary V-I curve corresponding to time t was constructed by repeating this for several clamping voltages. We note that momentary V-I curves can be computed for the general case of a pair of coupled neurons, and thus may have some advantage over phase plane methods.

	RESULTS

Abstract Introduction Methods Results Discussion References

Cells with different channel densities show various responses to current injection

Our neuron model could be approximately classified, depending on the calcium and leakage conductances, into four different types with qualitatively different response for steady current injection: a stable neuron, a spontaneous oscillator (SO), a conditional oscillator (CO), and a conditional bistable neuron. These four behavioral types are illustrated by the V time courses under current clamp in Fig. 2, insets. A stable neuron responds to current injection by settling to a unique and stable membrane potential, perhaps after a transient nonlinear response. The steady-state I-V relation is monotonic increasing. An SO is an autorhythmic neuron, requiring no input to oscillate. In general, it also oscillates for a range of applied currents, extending from some negative minimum to some positive maximum. A CO is driven to oscillate by injection of either strictly positive or strictly negative currents. Except for a very small range of parameters, our model cell is a CO only for negative current. We define a conditional bistable neuron as one that, in response to a range of currents, converges to one of two membrane potentials. The final membrane potential depends on the initial conditions, i.e., on the values of V and h before the current injection. Also, a brief current pulse can switch the membrane potential from one stable value to the other (not shown). Figure 2 shows the regions (separated by ) in the _T-g_L plane that correspond to the four neuron categories. The dotted line shows the conductance values at which the stable steady state changes from purely exponential decay to a damped oscillatory decay; that is, where the eigenvalues change from real to complex (see METHODS). Below this curve, our cell model has an oscillatory component (either sustained or damped). Thus some stable cells are damped oscillators. The divisions in Fig. 2 depend on the parameters (V_L and V_Ca) and on the gating kinetics for I_T. The techniques used to compute this behavioral segmentation are described in APPENDIX A.

View larger version (21K): [in this window] [in a new window]   FIG. 2. Different combinations of channel densities yield four types of cells. The T-gL plane is approximately divided into four different zones (separated by ), corresponding to model cells with different response properties to steady current injection. The zone labeled "stable" corresponds to cells whose membrane potential evolves to a stable steady state, unique for each Iapp value. Inset: voltage response of a stable cell (T = 0.4 mS/cm2, gL = 0.25 mS/cm2) to current pulses of 0.5, 0, and 0.5 µA/cm2 injected for 500 ms. Scale bars: 0.5 s, 5 mV. Cells within the "conditional bistable" zone converge to either 1 of 2 stable steady membrane potentials, depending on the initial conditions. This is exemplified by a cell with T = 0.4 mS/cm2 and gL = 0.06 mS/cm2 (inset). Current pulses of 0, 0.03, and 0.06 µA/cm2 are superimposed on a tonic current of 0.3 µA/cm2. In the case with Iapp = 0.06 µA/cm2, V does not return to its original level after the pulse terminates. Between these zones are cells that oscillate for some range of negative currents, the "conditional oscillator" zone. Inset: voltage time course of such a cell (T = 0.4 mS/cm2, gL = 0.1 mS/cm2), when injected with step currents of 0.18 (periodic activity), 0, and 0.18 µA/cm2. Cells in the "spontaneous oscillator" (SO) zone generate oscillations without any stimulus. Inset: membrane potential time course of such a cell (T = 0.4 mS/cm2, gL = 0.15 mS/cm2) with no current (periodic activity), as well as with Iapp = 0.5 and 0.5 µA/cm2, where the oscillations are abolished during the current pulse. Dotted line: curve where the eigenvalues change from real to complex (see METHODS). Below this curve, the eigenvalues are complex and the system has oscillatory components, either sustained or damped. See Appendix A for methods to determine zonal boundaries.

Note that the zones classifying the neurons by their response to current stimulus appear continuous on the _T-g_L plane, except for small parameter regions near zonal boundaries where mixed behaviors might occur. It is not essential for us to classify separately these mixed behaviors here (see APPENDIX A). When _T is very small, Eq. 14 implies that the condition for instability cannot hold, even if g_L is reduced proportionally. With larger values of _T, a neuron can change from a stable cell to an SO, then to a CO and finally to a conditional bistable cell by continuously decreasing g_L (or increasing _T).

In Fig. 3, A-D, we show the stimulus-response diagrams for the four types of cells mentioned above. Only g_L varies (decreasing) from A to D; the insets show the _T-g_L combination for each case. The diagrams represent the voltage only for steady (time-independent) or periodic states, both stable (solid) and unstable (dotted). The curve corresponding to the steady state thus constitutes a cell's steady-state I-V relation. A stable cell (A) has a unique stable steady-state V that increases monotonically with applied current. Its I-V relation is nearly linear except near the resting potential (I_app = 0), where it bends because of the existence of a steady-state ("window") calcium conductance. An SO (B) is characterized by having a stable limit cycle solution in a range including positive and negative currents. At any of these currents, the cell oscillates with voltages extending between the low and high values on the solid portion of the bubble-shaped curve. In the case of a CO (C), a stable limit cycle solution exists in a range of currents that is strictly negative or, as occurs in some small parameter range for our chosen I_T kinetics, strictly positive. For low enough g_L, the I-V relation develops an N-shaped character, and thus multiple stable steady states may coexist for some current levels. Figure 3D shows such a case, in which the cell is bistable. In a limited range of negative currents, two stable steady solutions coexist. Therefore, depending on initial conditions, the cell settles to either of two stable steady potentials. In the conditional bistable case, if a limit cycle exists for some range of applied current it is unstable in most cases. Such a cell cannot oscillate for any steady injected current. Although this empirical observation holds for our particular cell model, and may hold for some other models, one should not expect it (or other features shown in Fig. 2) to hold for behavioral state categorizations of excitable membrane systems in general.

View larger version (16K): [in this window] [in a new window]   FIG. 3. Response diagrams for current injection for the four types of neurons. The steady-state and periodic (repetitive firing) solutions, voltage amplitude vs. applied current, for four different cell types. Each inset shows the combination of calcium (T) and leak (gL) conductances (marked × on the T-gL plane copied here from Fig. 2) of that panel's cell. Solid and dotted/dashed lines: stable and unstable solutions, respectively. A: steady-state solution of a stable cell (T = 0.4 mS/cm2, gL = 0.25 mS/cm2). Resting potential of this cell(intersection with ordinate) is 61 mV. B: response states of an SO (T = 0.4 mS/cm2, gL = 0.17 mS/cm2). For Iapp less than 0.137 µA/cm2 and Iapp greater than 0.058 µA/cm2, a unique and stable fixed point solution exists. For 0.132 µA/cm2 < Iapp < 0.058 µA/cm2, a stable limit cycle solution coexists with an unstable fixed point solution. For example, with no current the cell oscillates between 54.3 and 60.3 mV (with a frequency of 5.4 Hz, not shown) with an unstable limit cycle (heavy dashes). C: current-voltage (I-V) curve of a conditional oscillator (T = 0.4 mS/cm2, gL = 0.11 mS/cm2) is N shaped. A stable limit cycle coexists with an unstable steady state for 0.284 µA/cm2 < Iapp < 0.114 µA/cm2. For other current levels there is a unique, stable steady state. Resting potential is 53.6 mV. D: case of a bistable cell. Channel densities are T = 0.4 mS/cm2, gL = 0.05 mS/cm2. For 0.434 µA/cm2 < Iapp < 0.235 µA/cm2, 2 stable steady states coexist. With no current, the fixed point (resting potential) is 48 mV. Note change in ordinate scale between A, B and C, D. LP1 and LP2: limit points 1 and 2.

Coupling two nonoscillating neurons may generate sustained oscillations

We begin our study of rhythmogenesis in olivary networks by considering in this paper the behavior of two coupled cells, identical in all but their channel densities (g_L and _T). For each of the two cells i,j = 1,2, the differential equation for the voltage is now

(15)

where g_coup is the coupling conductance. Intuitively, if two SOs are coupled electrically they will oscillate, phase-locked together, if their intrinsic frequencies are not very different. More interestingly, we show that cells that do not oscillate spontaneously when isolated may form sustained, low-amplitude oscillations when electrically coupled.

The most straightforward case to consider is a pair of cells that is strongly coupled. The limiting situation is for infinite coupling conductance (g_coup = ), and this is equivalent to a single neuron whose g_L and _T values are the average of the corresponding conductances in the two neurons. Clearly, only when these average conductances fall within the SO zone is the "average cell" oscillating spontaneously. If the average cell does oscillate, then one expects that the pair will also oscillate for some range of g_coup that extends from some finite value to infinity (and our analysis in APPENDIX B confirms this expectation for strong coupling). This observation, coupled with numerical simulations, gives us a large parameter range for generating oscillating cell pairs. Figure 4, A-E, shows the temporal behavior of five neuron pairs. For each pair, a few coupling conductances are examined (rows). The schematic above each column illustrates the type of neurons defined on the _T-g_L plane (see Fig. 2). For those pairs that are not oscillating at steady state, we inject a brief negative current (at ) to show how the system converges back to its steady, quiescent state. In Fig. 4, A-D, one of the two neurons is a stable cell (circle) and the other is a CO (square). In all four cases, both neurons are quiescent when isolated (g_coup = 0). However, transient hyperpolarization induces a short sequence of damped oscillations. When electrically coupled, the behaviors of the four pairs differ considerably. Figure 4A illustrates the concept stated above. The average of the two cells' channel densities (labeled with a cross) falls within the SO zone, and the pair forms sustained oscillations when g_coup is above some minimum value. These oscillations are low amplitude and in phase; the amplitude rises and the frequency decreases with increasing coupling conductance (not shown). For g_coup less than the minimum value, the cell pair shows damped oscillations following transient current perturbations. In Fig. 4B, the pair oscillates indefinitely for coupling conductances of 0.1 and 0.25 mS/cm², but not for a coupling conductance of 0.5 ms/cm². In this case, the amplitude rises and decreases with increasing coupling conductance, whereas the frequency decreases with increasing coupling conductance. Because this pair's average channel densities falls in the stable (ST) zone, we expect that with high coupling strengths this pair will be quiescent. The strongly coupled pair shows damped oscillatory behavior in response to transient hyperpolarization. Figure 4, C and D, demonstrates two cases that do not form sustained (only damped) oscillations at any of the coupling conductances shown, nor for any other coupling conductances (not shown).

View larger version (20K): [in this window] [in a new window]   FIG. 4. Time domain behavior of electrically coupled pairs of neurons. Shown are several selected examples of electrically coupled pairs of neurons. Each column (A-E) represents a specific pair of neurons; as in Fig. 3, the inset gives the corresponding cells' conductances (T, gL). Rows: different values of electrical coupling conductance gcoup, labeled at far right. Each pair of sweeps shows the two cells' voltage time courses. In pairs that were not spontaneously oscillating, a 20-ms, 0.1-µA/cm2 pulse was injected (at the time marked ). Squares and circles in the inset correspond to the top and bottom sweeps, respectively. Crosses in insets correspond to the "average" cell. A-D: 1 of the neurons is stable and the other is a conditional oscillator. In all four of these cases the two cells are quiescent when isolated (gcoup = 0 mS/cm2). Nevertheless, very different behaviors are observed when the cells are electrically coupled. In A, the pair forms sustained oscillations for electrical coupling of 0.25 and 0.5 mS/cm2. As the coupling increases, the amplitude increases, and the frequency decreases. Note that the average cell is an SO. B: pair oscillates indefinitely for electrical coupling of 0.1 and 0.25 mS/cm2, but not for higher coupling values. The average cell is a stable cell. C and D: two examples of pairs in which the oscillations damp out, no matter what coupling value is used. Note that in C the line connecting the two cells passes through the SO area, but the average cell is a conditional oscillator. D: line connecting the two cells does not cross the SO zone. E: case of a stable neuron electrically coupled with a bistable neuron. The bistable cell may settle to either of two membrane potentials in a negative range of applied currents. Oscillations are not sustained for coupling values of 0.0, 0.1, and 0.25 mS/cm2. Yet, with a coupling of 0.5 mS/cm2, the pair oscillates indefinitely. Note that the average cell of this pair is an SO.

Note that in Fig. 4C the average channel densities falls in the CO zone (just below the border of the SO zone). Note also that in Fig. 4D the line connecting the two cells on the _T-g_L plane does not cross the SO zone. Below we refer to these two observations in formulating a heuristic criterion to predict which pairs of neurons are likely to develop oscillations when electrically coupled. Figure 4E demonstrates a case where one of the two neurons was stable (circle) and the other was a conditional bistable cell (square). At steady state, the two cells are quiescent for coupling conductances of 0, 0.1, and 0.25 mS/cm², yet they generate a periodic sustained behavior with a coupling conductance of 0.5 mS/cm². The combination of average channel densities falls within the SO zone.

Except for the strong coupling case, it is generally difficult to predict analytically which coupled neuron pairs will spontaneously oscillate. Nevertheless, on the basis of our empirical experience we have developed the following "rule of thumb". We found that, in the majority of the cases where neither of the two isolated neurons is an SO, the pair generates sustained oscillations if the two following conditions hold.

1) The line connecting the two cells on the _T-g_L plane crosses the SO zone.

2) The combination of average channel densities falls inside the SO zone or in the ST zone. When the average is in the SO zone, the two cells oscillate for any coupling conductance above a critical value; when in the ST zone, the pair is oscillatory only for a finite range of coupling conductances.

[Of course, if either isolated cell is an SO, lying inside the SO zone, then the pair will oscillate if weakly coupled, that is, for some (perhaps small) range of g_coup, starting from g_coup = 0.]

To support our approximate rule, we systematically tested all possible pairs, with calcium and leak conductances assigned from the discrete range 0.05 + i * 0.05, i = {0,1, . . . ,9}. Thus the total number of pairs tested was 5,000. For each pair, we ran AUTO with the coupling conductance as the control parameter, ranging from 0 to 10 mS/cm². A pair was declared as oscillatory if one or two HB points were detected by AUTO. We found no combinations of non-SOs that could oscillate when their average is in the CO zone or in the "conditional bistable" zone. We found only two pairs (of the 5,000 tested) that could oscillate when the connecting line does not cross the SO zone, and in those two cases the connecting line passes very close to the zone's tip.

Next we study in detail the mechanism by which two nonspontaneously oscillating neurons generate oscillations when they are electrically coupled.

Mechanism for generating sustained oscillations in coupled neurons

Two different approaches are used to understand this gap-junction-mediated rhythmogenesis. In the first approach, we study extreme parameter regimes for which the full system of four ordinary differential equations can be approximated by analytical (asymptotic) methods to reduced systems of equations that can be more easily interpreted and solved. In the second approach, we study stability empirically via momentary I-V relations, by simulating the full system with one of the two cells voltage clamped.

For the case of strong coupling, we use perturbation theory to develop an approximate description for a pair's behavior. To implement this method it is useful to define, for comparison purposes, a reference intrinsic conductance g_ref; for example, we could use g_ref = g_L. In the approximation scheme the ratio g_ref/g_coup is treated as a small parameter  (see APPENDIX B for details, and see Nayfeh 1973 for a general description of perturbation methods). We assume that each cell's variables can be represented as a perturbation series in powers of . The lowest-order term for the membrane potentials (corresponding to the case of infinite coupling) is V_av, the membrane potential of the average cell. Thus

(16)

The subscript j here distinguishes the two cells, j = 1,2. The second subscript denotes the coefficients for the successive terms in the series approximation. These coefficients are of order one, that is, independent of . With the use of a corresponding series for the h variables, one then substitutes these representations into the four differential equations. Adding the two current balance equations and setting  to zero yields the equation for V_av

(17)

where the terms I_ion,j contain the different values of _T and g_L for the two cells; the factor 2 multiplying C_m comes from having added the two equations and yields the expected average. The equations for h_av in the two cells are identical to lowest order. Next we show that there is a significant coupling current that flows between two nonidentical cells, even when they are strongly coupled (a reader who wishes to bypass the following derivation may advance to Eq. 21 and the succeeding text).

The coupling current that flows through the gap junctions from cell 2 to cell 1 is defined as

(18)

Inserting the series approximations into Eq. 18 and noting that the common terms V_av cancel, we find that to lowest order

(19)

In general, _coup is not small but of order one when the cells are nonidentical. Going back and subtracting the current balance equations for the two cells, using the series representations and keeping only the lowest-order terms, we get

(20)

where the voltage argument of I_ion,j is V_av. This equation describes approximately the dynamics of I_coup. The leading factor C_m/g_coup is the time constant for the relaxation of _coup, and it is very brief compared with the reference time constant C_m/g_ref.

From this equation, we see that I_coup equilibrates very rapidly to approximately half the difference (I_ion,1  I_ion,2). Then, after a transient phase (with duration of order _h) during which the gating variables equilibrate to h_av, we have

(21)

This result explicitly shows that if the cells are nonidentical, the coupling current is of size comparable with the ionic currents. For strongly coupled cells to achieve near-equipotentiality, I_coup acts as an equalizing current to offset the difference between their ionic currents, because of the different _T and g_L values. The result also shows that, generally, I_coup is time varying.

Some of these points are illustrated in Figs. 5 and 6 for a case similar to that in Fig. 4A. We start by describing the case of infinite coupling (Fig. 5A). The network composed of a stable cell (cell 1) infinitely coupled to a CO (cell 2) is equivalent in its dynamic behavior to a cell with the average channel densities. Thus this dynamic behavior can be evaluated by examining the nullclines of the average cell. The V nullclines of cell 1 (· · ·), cell 2 (- - -), and the average cell () are plotted in the V-h plane, along with the h nullcline (), which is the same for each of the three cases. The positive slopes at which the V nullclines of cell 1 and cell 2 intersect with the h nullcline show that these two cells are quiescent at rest, i.e., with no injection of current (see METHODS). However, the intersection of the V and h nullclines of the average cell occurs at a negative slope of the V nullcline (arrow). A simple calculation shows that with the parameters chosen for this case, Eq. 14 is satisfied. Thus the fixed point of the average cell is unstable; the average cell has a stable limit cycle solution, predicting that for strong coupling this pair should oscillate spontaneously. Such a case can be seen in Fig. 5B, where the two cells are electrically coupled with g_coup = 0.5 mS/cm².The limit cycle trajectories of cell 1 (· · ·) and cell 2(- - -), which were computed by simulation of the full system, are superimposed on the nullclines of the average system (). These periodic trajectories are quite close to the periodic trajectory computed for the average equation (). As must be, this limit cycle trajectory crosses the average cell's V and h nullclines horizontally and vertically, respectively. We note here that the amplitude of the h oscillation, as for the V oscillation, is quite small.

View larger version (15K): [in this window] [in a new window]   FIG. 6. Comparison of actual coupling current with approximated coupling current in the limit of infinite coupling. Cell 1 and cell 2 are the same as in Figs. 4A and 5. A: actual coupling current (Icoup, computed with Eq. 18) when the two cells are coupled with gcoup = 0.5 mS/cm2, and the coupling current predicted from the series approximation in the limit of strong coupling (coup, computed with Eq. 21, dotted line). B: stable solutions of the coupling current are plotted as functions of gcoup (). For values <0.13 mS/cm2, the stable solution is a fixed point. When gcoup  0.13 mS/cm2, the stable solution oscillates between the two values shown on the plot. Dotted horizontal lines: maximal and minimal values of the coupling current predicted by the series approximation in the limit of strong coupling.

Figure 5C shows the time courses of the ionic currents (I_ion,1, · · ·; I_ion,2, - - -) and the coupling currents (I_coup and I_coup, ), when cell 1 and cell 2 are electrically coupled with g_coup = 0.5 mS/cm².

From the definition given in Eq. 18, I_coup is the current flowing from cell 2 to cell 1. For each cell, the coupling current's time course is seen as comparable in amplitude with the ionic current. As our perturbation results show, this statement remains true even for very strong coupling. Curiously, in this particular case, all four currents oscillate as expected but they do not change sign: I_ion,1 and I_coup are strictly depolarizing, whereas the opposite is true for I_ion,2 and I_coup. However, the net membrane current, ionic plus coupling, does change sign: it periodically varies from positive (when, from Eq. 15, the voltage decreases) to negative (when the voltage increases). This can be seen by comparing, for example, I_ion,1 with I_coup. Whenever I_ion,1 > I_coup, the difference I_ion,1  I_coup is positive and dV₁/dt < 0; whenI_ion,1 < I_coup, dV₁/dt > 0.

Figure 6 illustrates the accuracy of the strong coupling approximation for I_coup derived in Eq. 21. In Fig. 6A, the time courses for the actual and approximate coupling current are compared for the case of Fig. 5. Even though the coupling conductance is not very large (g_coup = 0.5 mS/cm²), the period and minimum amplitude are quite accurate; the maximum amplitude differs by only ~20%. In Fig. 6B one sees that the accuracy improves as g_coup increases. Note, to lowest order, the approximate coupling current I_coup does not depend on g_coup; thus the approximate amplitude measures appear here as the horizontal dotted lines. The amplitude of I_coup asymptotically approaches that of the approximate coupling current _coup as g_coup increases, in a hyperbolic fashion as expected. Also in Fig. 6B we see that when g_coup < 0.13 mS/cm², the cell pair no longer oscillates. The emergence of the oscillation is an HB at the critical coupling conductance; for g_coup > 0.13 mS/cm², the system's fixed point (not shown) is unstable. Figure 6B shows explicitly the behavior for a cell pair for which the average falls in the spontaneously oscillating zone, i.e., an oscillation for g_coup above some critical value. Finally we note that, as can often happen with perturbation theory, the approximation can give accurate results even when  is not very small. In the above example, g_coup is only ~2-5 times larger than the ionic conductancesnot tens of times larger.

Another extreme case that can be used to simplify the system of differential equations is one where the coupling is moderate, but one of the two cells (for example, cell 1) has a very large g_L, compared with _T and g_coup. Thus cell 1 is a strongly stable cell. dV₁/dt is governed mainly by the difference between V₁ and V_L, and V₁ relaxes quickly to make this difference small, i.e., to achieve V₁  V_L. The other cell (cell 2) is a CO or a conditional bistable cell. The differential equation for its voltage reduces to

(22)

or

(23)

where g_L = g_L + g_coup. In other words, the effect of coupling cell 1 to cell 2 is approximately equivalent to a mere increase of leak conductance for cell 2.

From Fig. 2 it can be seen that if the "new" effective channel densities of cell 2 (_T and g_L) are within the SO zone, the system spontaneously oscillates. This implies that there is a range, a finite range, of g_coup for which this particular pair of cells forms sustained oscillations. This observation is consistent with the rule of thumb that was described in the previous section, because 1) if cell 1 is a "strong" stable cell, the average of cell 1 and cell 2 falls within the ST zone and 2) the line linking the two cells crosses the SO zone (in the strong g_L limit, the line is vertical and thus the oscillation occurs only if the _T value for cell 2 exceeds the minimum _T value in the SO zone).

Our rule of thumb can also be supported by considering a complementary limit. Suppose one of the two cells (for example, cell 1) has a very large _T, compared with g_L and g_coup. This cell is far to the right beyond the SO region. Therefore we expect V₁  V_Ca. Thus for cell 2 we have approximately

(24)

The last term in this equation can be rewritten as

(25)

Now inserting this into the current balance equation we get

(26)

where g_L = g_L + g_coup and I = g_coup (V_Ca  V_L). In this case, coupling these two cells is approximately equivalent to increasing the leakage conductance of cell 2 and to injecting a steady depolarizing current into cell 2. We conclude from this that if the isolated cell 2 was of the stable type, this pair will not oscillate, because the increased leakage conductance pushes it farther into the ST regime and depolarizing current (I here) cannot make a stable cell oscillate (by definition). This pair, whose average conductances fall below the SO zone, will not oscillate for any range of g_coupas our empirical rule of thumb states.

Consistent with our parenthetical statement succeeding our rule of thumb, if cell 2 was spontaneously oscillating and coupled weakly with cell 1 then the pair would oscillate. According to the above demonstration, if g_coup is small the conductances _T and g_L would still lie in SO zone and a small depolarizing current would not destroy the stable oscillation. On the other hand, the oscillation will not persist for nonweak coupling. In such a case, one of two silencing events will occur: either g_L increases to "move cell 2" into the stable regime or I increases to steadily depolarize cell 2 while it is still in SO.

One might argue intuitively about the mechanism for why a pair of cells (one stable and the other a CO or a conditional bistable neuron) could oscillate when coupled, as follows. The stable neuron, cell 1, with a more negative resting potential, could hyperpolarize cell 2, bringing it into the oscillating regime. This argument could be used for the extreme case of very large leak conductance in cell 1, as described in the previous paragraph. But in general the argument is incomplete. It ignores the converse behavior, that coupling would tend to depolarize cell 1, and contradicts the fact that cell 1 does not oscillate for steady depolarization. The argument is oversimplified; it can be sharpened on the basis of our perturbation results. One should not think of the coupling current as steady or time averaged. It is time varying, and in such a way that it transiently depolarizes cell 1 and transiently hyperpolarizes cell 2 to obtain the regenerative response of the synchronized oscillation.

In the case of nonstrong coupling and moderate values of the leak conductance, we have no analytic treatment for rhythmogenesis. We cannot compute the eigenvalues of the system in general, and phase plane diagrams cannot be drawn. However, one can gain some insights by computing momentary V-I curves for this case (see METHODS). In Fig. 7, momentary V-I curves at 3, 10, and 300 ms are plotted for a stable cell (cell 1, Fig. 7A), a CO (cell 2, Fig. 7B), an average cell (coupling of cell 1 and cell 2 with an infinite coupling, Fig. 7C), and cell 1 coupled to cell 2 with a finite coupling (Fig. 7D). In the case of a single cell (Fig. 7, A-C), the ordinate represents the ionic current. In the case of a pair of coupled cells (D), it represents the ionic plus coupling current. In this latter case, the clamped cell is cell 2. In the case of cells 1 or 2 (Fig. 7, A or B), the fixed point occurs on positive slopes of the momentary V-I curve. Thus a small depolarizing deviation from the fixed point results in a positive (outward) ionic current, which would hyperpolarize the membrane potential back toward the fixed point, in a nonclamped situation. Likewise, a small hyperpolarizing displacement leads to a negative (inward) current that depolarizes the membrane back toward the fixed point. In the case of the average cell (Fig. 7C), or the coupled pair of cells (Fig. 7D), there exists a period of time in which the momentary V-I curve has a negative slope at the fixed point (Fig. 7C, 3 and 10 ms; Fig. 7D, 10 ms). Thus, during this brief time window, a small depolarization generates an inward current (and a small hyperpolarization generates an outward current). If this current is sufficiently large and long lasting to produce a substantial voltage change, the membrane potential (in the nonclamped situation) will be driven away from the fixed point. A large inactivation time constant (_h) ensures that the destabilizing current will not be short lived. In such a case, the resting potential of cell 2 is not stable. Because of the electrical coupling, it entrains cell 1 to oscillate with it.

View larger version (14K): [in this window] [in a new window]   FIG. 7. Analysis of stability with momentary voltage-current (V-I) curves. Cells 1 and 2 are a stable cell and a conditional oscillator, respectively, as defined in Fig. 5. The momentary V-I curves for 3, 10, and 300 ms are shown for cell 1 isolated (A), cell 2 isolated (B), the average cell (C), and the pair, cell 2 and cell 1, coupled with gcoup = 0.1 mS/cm2 (D). In the last case, cell 2 is voltage clamped and the V-I curve shown is for cell 2. In the cases of cell 1 and cell 2, all 3 momentary I-V curves intersect at the fixed point with a positive slope. The fixed point voltages are 59.8 and 52.9 mV, respectively. In the case of the average cell, the momentary I-V curves for 3 and 10 ms intersect at the fixed point with a negative slope, as also occurs for the momentary I-V curve for 10 ms in the case of the coupled pair. In these 2 cases, the fixed point voltages are 56.6 and 56.1 mV, respectively.

	DISCUSSION

Abstract Introduction Methods Results Discussion References

Understanding the relative contributions of coupling and intrinsic properties to rhythmogenesis in networks of excitable cells requires combined modeling and experimental efforts. Some insights into how network oscillations are established and coordinated are being provided by case studies of a variety of systems, such as the pancreatic -cells (Meda et al. 1984), the pyloric network of the stomatogastric ganglion in crustacea (Hooper and Marder 1987), the central pattern generator of a swimming lamprey (Grillner and Matsushima 1991), the locus coerulus in neonatal rats (Christie et al. 1989), the IO complex (Llinás and Yarom 1986a), etc. Theoretical studies show that electrical coupling can modulate the frequency of the oscillations, either reducing or increasing it (Abbott et al. 1991; Kepler et al. 1990; Meunier 1992). Other studies show, counterintuitively, that weak coupling can produce antiphase synchronization with oscillatory, or even nonoscillatory, neurons (Sherman and Rinzel 1992), or that mutual inhibition can lead to synchrony when the postsynaptic conductance decays slowly (Wang and Rinzel 1992). Coupling of nonlinear elements yields complex dynamics, and there are many gaps in our levels of understanding. The complexity and richness of behaviors are sources of flexibility and plasticity, properties that have major importance in nervous systems.

The present work was inspired by observed behaviors in the IO complex, which is composed of neurons that oscillate in a subthreshold range (Benardo and Foster 1986; Llinás and Yarom 1986a) and that are extensively coupled electrically (de Zeeuw et al. 1990; Sotelo et al. 1974). The electrical coupling between IO neurons is not only important for synchronization (Llinás and Sasaki 1989; Sasaki et al. 1989), but probably plays a key role in the generation of the sustained oscillatory behavior of the neurons (Bleasel and Pettigrew 1992, 1994; Yarom 1989, 1991). Several lines of evidence suggest that the IO neurons, at least most of them, are not intrinsic oscillators. Their dynamic behavior is better described as that of damped oscillators. Electrical coupling per se cannot mediate sustained synchronous oscillations among such units, if they are completely identical. However, IO cells exhibit a significant variability in input resistance, resting potential, membrane time constant, magnitude of calcium current, and resonant frequency (Lampl and Yarom 1997; Manor 1995). Heterogeneity in ion channel densities may contribute to these variabilities. We hereby suggest that this heterogeneity may have functional importance for the generation of spontaneous oscillations. Heterogeneity in channel densities among classes of neurons is frequently observed in other preparations (Benitah et al. 1993; Christ et al. 1993; Hájos and Greenfield 1993). Theoretical models of networks of excitable cells showed that such heterogeneities may play a functional role (Smolen et al. 1993). The extensive electrical coupling and the large electrical variability observed in the IO inspire the question: can rhythmogenesis be expected from electrical coupling of neurons with different channel densities, even if none, or just a few, of the individual cells are SOs? Before answering this question, we describe the behaviors that should be expected from an isolated IO cell as a function of its channel densities.

We used a biophysically based minimal cell model with only two conductances (_T, g_L) to map out these behaviors. This simplification of the IO neuron was motivated by the experimental findings that neither potassium nor sodium nor high-threshold calcium channels are required for the existence of the STO (Benardo and Foster 1986; Lampl 1994; Llinás and Yarom 1986a). We found that differences in channel densities yield a spectrum of behaviors that can be categorized into several distinct types of response to current injection (Fig. 2). We characterized four such major behaviors: stable cells, SOs, COs (which require current injection to oscillate), and conditional bistable cells. The channel densities (and gating properties of the T-type calcium current) interact to significantly affect the neuron's "excitability." First, the leak conductance is the major determinant of the input resistance. In this sense, it acts as a stabilizing parameter: the larger it is, the smaller and "more passive" are the cell's responses to input or intrinsic currents. Second, the relative values of the leak conductance and the maximal calcium conductance determine the resting membrane potential. In our case, where the calcium current has some window conductance, the value of the membrane potential relative to the window conductance region determines the excitability of the cell in a dramatic way. For example, instability may occur only when the resting potential is in a range where the short-time transient slope conductance is negative. Such a condition may exist only if the voltage-dependent conductance is partially open at rest and not completely inactivated, i.e., if the resting membrane potential is within the window conductance range. Third, changes in the maximal calcium conductance modify the relative importance of the window conductance. Thus it acts as a destabilizing parameter (in the range of voltages of the window conductance). Note that the extracted kinetics of the low-threshold calcium current from IO neurons (Fig. 1A) suggest that a window conductance does exist. Moreover, the dependence of resonance (existence of a range of input frequencies, where the voltage response is maximal) on a calcium window conductance (Hutcheon et al. 1994), and the demonstration of resonance in the IO neurons (Lampl and Yarom 1996), further support the existence of a calcium window conductance in the IO neurons. Fourth, the time scale of I_T's inactivation gating is absolutely critical in determining stability of the rest state. Although the window conductance is essential, it is a steady-state property and its existence does not guarantee an oscillation. The inactivation must be slow enough relative to the destabilizing time scale of the negative resistance.

Our model, which is partly motivated by observed parameter variability in slice preparations, suggests that there are at least four different behavioral types of IO neurons. Not all these different types have been identified in experiments. Rather, with the exception of spontaneously oscillating neurons found in 10% of the slices (Llinás and Yarom 1986a), most intracellular recordings from slices showed stable quiescent behavior with damped oscillatory transients. These experimental findings, however, do not preclude the possibility that the other types exist in the population of IO neurons. One expects that behavioral properties of the individual neurons are "blurred" because of electrical coupling with their neighbors, as we have found in a large network model (unpublished results). Our classification as described in Fig. 2 may only be feasibly attempted experimentally after a highly specific blocking agent for gap junctions becomes available. The experimental finding that most stable IO neurons generate damped oscillations after brief intracellular current injection is supported by our simulations with the two-neuron model (cf. Fig. 4) and with large networks (unpublished data).

From the above characterizations we learned that, for most values of the channel densities in Fig. 2, the individual model neurons are not spontaneously oscillating. Our results show, however, that coupling two quiescent cells can lead to a pair showing spontaneous oscillations. This occurs, for proper coupling ranges, if the individual neurons have different channel densities whose values straddle the SO zone in Fig. 2 and whose average values lie in or above this region. The electrotonic coupling (if not weak) creates a two-cell network that displays a behavioral pattern different from that of the individual cells. The source of this new, emergent behavior lies in the nonlinear effect of the equilibrium potential of each separate cell on the calcium conductance. For example, consider the case of two neurons, one with a resting potential that is more hyperpolarized and the other that is more depolarized compared with the voltage range where the window conductance exists. Both neurons are stable, but their electrical coupling can shift their voltages into the window conductance regime, where their calcium conductances are open and may destabilize the pair. Another surprising example is the behavior of a pair of neurons composed of a conditional bistable neuron coupled to a stable neuron. Neither of the two neurons is capable of generating a periodic behavior with any current injection. Yet, electronic coupling with the appropriate coupling strength may yield sustained oscillations in the two-neuron network (Fig. 4E). Our analytic perturbation theory for the case of strong coupling has shown us that the electronic coupling current is the mediating mechanism for new behaviors when the cells are not identical. Even in the limit of infinite coupling, this current is sizable although the cell voltages are nearly identical. The coupling current makes up the difference in the two cells' ionic currents caused by their different channel densities.

The idea that electrical coupling can contribute to rhythmogenesis in a heterogeneous population was previously suggested for bursting pancreatic -cells. Smolen et al. (1993) postulated that bursting oscillations could arise by parameter averaging among electrically coupled quiescent and continuously spiking cells. Our analytic result (cf. APPENDIX B) shows that averaging is the rhythmogenic mechanism in the case of strong coupling. This result is general. Any pair of nonidentical cells that can be modeled as in Eq. 21 (even if there are more ionic currents and gating variables) will oscillate for sufficiently large g_coup if the average cell oscillates spontaneously. The result is restricted in that we consider heterogeneity only in the I_ion expressions, thus allowing differences in the channel densities and reversal potentials but not in gating or other kinetics. In contrast to the study by Smolen et al. (1993), where variability was described in terms of an individual cell's spontaneous activity, our categorization is based on how intrinsic parameters change the cell's stimulus-response properties for current injection. That is, Fig. 2 distinguishes behavior for ranges of three, not just two, parametersthe two channel densities and the stimulating current. Although in both studies oscillations could sometimes arise even if the mean parameters fell outside the SO zone, our rule of thumb describes the conditions under which this could happen (for our model). Whether this aspect generalizes, we do not yet know.

From our work showing that heterogeneity in channel densities may generate new behaviors via coupling of the neurons and from the study by Smolen et al. (1993), we expect that other types of variabilities might yield similar results. Differences in kinetics of the voltage-dependent conductance (Berlind 1993) or in local concentrations of extracellular potassium (Guckenheimer and Labouriau 1993) are likely to produce different neural behaviors. These differences, via coupling, can also support network behaviors that are not expected from the response of its individual elements. Differences in cell sizes may also produce new behaviors via coupling.

Our study of heterogeneity and electrical coupling with the two-neuron model paves the way for a more detailed investigation and model of the olivary nucleus. The model we present obviously has some limitations and should not be applied to address certain issues. For example, because only two neurons are coupled together, intracellular injection of current to one of the cells can generate STOs, change their frequency, or abolish them. In experiments, the STOs are insensitive to current injection in a single cell. The basic components of our model, when incorporated into a large network, lead to behaviors that do match with such experimental observations. For example, the frequency of STOs in different slices is variable, between 3 and 10 Hz; this frequency range is seen in simulations with multicellular network models, where a moderate number of heterogeneous cells were used (unpublished results). Here, we have developed the framework for this, addressing the next set of questions concerning multicellular behavior in the IO nucleus.

From a generalist's point of view, our simplified model of the IO provides additional insights into the emergence of electrical behavior of isolated cells and electrically coupled neurons. We have shown that heterogeneity in channel density, integrated through electrical coupling, can result in a rich repertoire of electrical activity in the subthreshold regime. By virtue of electrical coupling, oscillatory behaviors can emerge from cells that are silent otherwise. Indeed, being alone can be very different from being coupled.

	ACKNOWLEDGEMENTS

  The authors thank V. Booth, B. Ermentrout, D. Hansel, D. Golomb, I. Lampl, and A. Sherman for helpful discussions and fruitful remarks.

  This work was supported by U.S. Office of Naval Research Grant N00014-19-J-1350 and by the United States-Israel Binational Science Foundation (Y. Yarom).

	APPENDIX A: SEGMENTATION OF THE T-gL SPACE INTO BEHAVIORAL "ZONES"

In this section we briefly describe the methods that we used to identify the regions in the _T-g_L plane that correspond to stable cells, SOs, COs, and conditional bistable cells.

First we explain how the SO zone (the wedge-shaped region in Fig. 2) is found. With g_L fixed, with the use of AUTO we compute the steady-state and periodic solutions as functions of the parameter _T; the results are presented in the bifurcation diagram of Fig. 8A. We start with _T = 0, where the steady state is = V_L, = h(V_L). AUTO computes the steady-state solution(s) by varying _T continuously over a preset range. It determines the solution's stability (stable plotted as solid line, unstable as dotted line) by inspecting the real parts of the eigenvalues. The transition points (between stable and unstable) are identified as HB₁ and HB₂. These are points from which periodic solutions of small amplitude emerge. In a second pass, AUTO can start from either of these points and compute the limit cycle solution as _T is varied until the other HB point is reached. Next, to get the "wedge," we use AUTO's two-parameter continuation feature. This allows us to track the two HB points as both _T and g_L are allowed to change (Fig. 2B). A cell with (_T, g_L) lying in the wedge has an unstable steady state and oscillates spontaneously. We remark that in some small parameter ranges we found that an HB might be subcritical, i.e., the emergent branch of limit cycles locally enters the regime where the steady state is stable before bending back into and traversing the wedge.

View larger version (10K): [in this window] [in a new window]   FIG. 8. Computation of the SO zone. A: voltage amplitude of the steady response shown as a function of the calcium conductance (T) (obtained by one-parameter continuation with AUTO). Leak conductance is gL = 0.3 mS/cm2. Solid and dotted lines: stable and unstable solutions, respectively. At T = 0.637 mS/cm2 and T = 0.936 mS/cm2, 2 HBs occur (labeled HB1 and HB2). B: 2 HBs were traced (with an AUTO 2-parameter continuation as both T and gL are allowed to change) and shown to be 2 branches of a continuous curve. At each gL value, the T values on the HB1 and HB2 branches define the calcium conductance range at which the cell is oscillating spontaneously. Note that with gL < 0.096 mS/cm2 or T < 0.237 mS/cm2, a cell cannot oscillate spontaneously.

The next stage is to define the zones for a stable cell, a CO, and a bistable cell. For a specific leak conductance (g_L), we start with a cell for which the ratio _T/g_L is large. Such a cell will have an N-shaped steady-state I-V relation, and perhaps be bistable, having two stable steady states for a range of currents. With AUTO, we compute the one-parameter bifurcation diagram: steady-state and periodic solutions as a function of I_app (Fig. 9A). Note that in this case the HBs are subcritical and the periodic solutions are unstable. In addition to the HB point(s) (1 or 2), AUTO identifies the two limit points LP₁ and LP₂ (where the steady-state solution plot has vertical slope). Now, with two-parameter continuation in AUTO, we trace these points as both I_app and _T are allowed to change (Fig. 9B). The curve of HB points is plotted as a solid line; the curve of limit points is shown as a dashed line. Note that, by definition, the voltages at which HB₁ and LP₁ occur are lower than those at which HB₂ and LP₂ occur, respectively.

View larger version (9K): [in this window] [in a new window]   FIG. 9. Method of classifying cell response types for a given gL. The gL value in this example (A and B) is 0.3 mS/cm2. A: steady-state solution(s) of the voltage shown as a function of Iapp (AUTO computation). In this case, T = 2 mS/cm2. Thick dashed lines: maxima and minima of unstable limit cycle solutions. Thin solid and dotted lines: stable and unstable fixed point solutions. By definition, the voltage at HB1 is lower than that at HB2. The same applies for LP1 and LP2. For 1.491 µA/cm2  Iapp  1.286 µA/cm2 (between HB1 and HB2), two stable fixed point solutions coexist. B: two HBs and the two LPs were traced with an AUTO continuation of two parameters (both Iapp and T are allowed to change). g0 is the T value at which the two HBs coalesce (0.636 mS/cm2). g1 andg2 are the two T values at which the HB curve crosses Iapp = 0 (0.6378 and 0.936 mS/cm2, respectively). g3 is the T value at which the two HBs intersect, i.e., the lowest T for which the current at HB1 is more positive than the current at HB2 (1.811 mS/cm2). Cells with T < g0 are stable (arrow at right). The range of values g0 < T < g1 correspond to cells that oscillate with injection of positive current only. In this figure, this range is very narrow and practically undetectable. When g1 < T < g2, the cell is oscillating with no injection of current (SO). Cells with g2 < T < g3 are oscillating in a negative range of currents and are thus called conditional oscillators. Finally, cells with T > g3 are "bistable."

We define four critical values of _T: the value where the two HB points coalesce g₀; the two values g₁, g₂ where the HB curve crosses the axis I_app = 0; and where the HB₂ and HB₁ curves coalesce, g₃. Note, the values g₁, g₂ will not exist if the HB loop lies strictly in the region I_app < 0, i.e., if the SO zone in Fig. 2 lies above the horizontal level for the specified g_L value. With the use of these critical values we classify cells with this g_L value in four different categories.

1) By the definition of g₀, at any _T smaller than g₀ there are no bifurcation points at any injected current. Thus the steady-state solution is unique and always stable, for any I_app. Thus the range _T < g₀ defines the stable neurons.

2) When g_1  T  g₂, the fixed point is unstable and the cell oscillates with no current; this defines the SOs.

In general, there will be a small regime where g_0  T < g₁, which corresponds to a sliver in the _T-g_L plane between the ST zone and the SO wedge.

3) When g₂ < _T < g₃ (or, in case the HB loop does not cross the _T axis, when g₀ < _T < g₃), bistability cannot occur because I(V_HB2) > I(V_HB1). g₂ marks the lowest value of _T for which both HB points occur for negative current. Thus these cells that oscillate only in a negative range of currents are called COs.

4) The merging of HB₁ and HB₂ at g₃ marks the lowest _T value for which I(V_HB2)  I(V_HB1). Thus, for _T  g₃, there is a range of negative currents at which two stable steady-state solutions exist. These are the conditional bistable cells.

These previous steps are repeated for different values of g_L. The g₀ values found for different g_Ls define the border of the ST regime. Note that when g₀ is plotted as a function of g_L, it defines the boundary between the ST zone and the CO zone in Fig. 2. It can also be seen that g₁ and g₂ as functions of g_L define the left and right boundaries, respectively, of the SO wedge in Fig. 2. The g₃ values found for different g_Ls mark the border segregating bistable cells from other cells in the _T-g_L space. Between g₀ and g₃, for a given g_L, are located the oscillators, spontaneous and conditional. Their distinction arises after we plot g₁ and g₂.

There might be some small overlaps between regions (or cases that are mixed) whose consideration would require more precise and complicated mathematical considerations than are worthwhile for our purposes. For example, a mixed case could be imagined as a distortion of Fig. 9A in which the unstable periodic branch from HB₂ might start as shown and then stabilize by bending back to more negative I_app; this cell would have for some small range of I_app a stable steady state (with V around 70 mV) coexisting with a stable limit cycle (with average V around 50 mV).

	APPENDIX B: APPROXIMATION OF THE COUPLING CURRENT IN THE CASE OF STRONG COUPLING WITH PERTURBATION ANALYSIS TECHNIQUES

We are interested in obtaining an approximation for the coupling current when the coupling conductance is strong (or, equivalently, when the coupling resistance is small). To do that, we ask how the solution is determined when the coupling resistance (the parameter) is slightly altered from zero. This type of analysis can be carried out with perturbation techniques. According to these techniques, the solution can be approximated by the first few terms of an asymptotic expansion, which is done in terms of the small parameter. Usually, the use of the first two terms gives an acceptable approximation.

In a proper perturbation treatment, the variables are first written in dimensionless form. This allows one to expose the small parameter that is to be used in the perturbation series. For the nondimensionalization, g_L, for example, could be defined as the reference conductance, g_ref. Then, the current balance equation (Eq. 15) is divided by g_ref and I_ion is redefined as The reference time constant is defined as The small parameter is then Similarly, t is scaled relative to _ref, and is now For each of the two cells j,i = 1,2, the current balance equation is now

(27)

For convenience, the gating equations (Eq. 5) are written as

(28)

To obtain an approximate representation when  is small (i.e., when the coupling conductance is strong), the voltages of the two cells are written as a perturbation series in powers of

(29)

where the (time-varying) coefficients _j,0, _j,1,··· remain to be found.

In problems that involve nonlinear oscillations, truncated expansions to the above form that do not account for the effect of  on the period are valid only for a finite time; this is because resonant behavior eventually leads to growth of the coefficient terms. In the Lindstedt-Poincaré method, these resonant behaviors are controlled by allowing the unknown frequency of the oscillations, ₁, to also depend on  (Nayfeh 1973). Thus this frequency is also expanded in series form

(30)

where ₀, ₁, etc. are to be determined. To introduce  into the problem, t is scaled by  so we obtain a new time variable The differential equations of our problem now take the form

(31)

and

(32)

For convenience, we define the following expressions Now, the system of equations defined in Eq. 31 and 32 is rewritten by adding and subtracting the differential equations for the voltages, and doing the same for the gating variables. After dividing by 2, we get

(33)

(34)

(35)

(36)

From Eq. 29

(37)

Similarly, , , and are expressed as series representations. These series representations, together with Eq. 30, are inserted into Eq. 33-36. For example, for Eq. 33

(38)

The next step is to find the lowest-order terms (₀, ₀, etc.). By setting  = 0 in Eq. 34, we get

(39)

This means that at the lowest order, the two voltages are identical. We denote this coequal voltage as _av(s)

(40)

We note that, by definition, ₀ = _av. By setting  = 0 in Eq. 35 and 36, we get

(41)

and

(42)

Using the definition of H(V,h) {= _ref[h(V)  h]/_h} in Eq. 42, evaluating, and then converting back to physical time, we see that Eq. 42 reduces to

(43)

Thus ₀ decays to zero with time constant _h. That is, to lowest order, the gating variables, if not equal at t = 0, equalize on the time scale _h. This allows us to define their equilibrated value It follows that, after this equilibration phase, ₀ = h_av, where according to Eq. 41 h_av satisfies

(44)

Finally, setting  = 0 in Eq. 38 and noting that ₀ = _av, we get (again, after a few time constants _h)

(45)

One could proceed to get the next-order terms, but we will not carry this to completion. However, it is valuable to consider the first step, because we get for "free," without having to solve additional equations, a useful approximation for I_coup. Notice that, because ₀ = 0,  = ₁(s) + ²₂(s) + ··· . We insert this expression into Eq. 34 and divide by . Setting  = 0 gives

(46)

By definition, I_coup = 2/ = 2(₁ + ₂ + ···). Thus with = 0

In other words, the coupling current needed to equalize the voltages of the two cells is simply proportional to the difference in their ionic currents.

The other unknown first-order terms (₁, ₁, and ₁) can be obtained by expanding the nonlinear terms in powers of . Note that in solving for these terms it is absolutely necessary that we have introduced the frequency into the perturbation series, to prevent these terms from growing unbounded in time.

	FOOTNOTES

  Address for reprint requests: Y. Manor, Volen Center for Complex Systems, Brandeis University, 415 South St., Waltham, MA 02254.

  Received 19 April 1996; accepted in final form 2 January 1997.

	REFERENCES

Abstract Introduction Methods Results Discussion References

• ABBOTT, L. F., MARDER, E., HOOPER, S. L. Oscillating networks: control of burst duration by electrically coupled neurons. Neural Comput. 3: 487-497, 1991.
• ALBUS, J. S. A theory of cerebellar function. Math. Biosci. 10: 25-61, 1971.
• BENARDO, L. S., FOSTER, R. E. Oscillatory behavior in inferior olive neurons: mechanism, modulation, cell aggregates. Brain Res. Bull. 17: 773-784, 1986.[Medline]
• BENITAH, J.-P., GOMEZ, A. M., BAILLY, P., DA-PONTE, J.-P., BERSON, G., DELGADO, C., LORENTE, P. Heterogeneity of the early outward current in ventricular cells isolated from normal and hypertrophied rat hearts. J. Physiol. Lond. 469: 111-138, 1993.[Abstract/Free Full Text]
• BERLIND, A. Heterogeneity of motoneuron driver potential properties along the anterior-posterior axis of the lobster cardiac ganglion. Brain Res. 609: 51-58, 1993.[Medline]
• BLEASEL, A. F., PETTIGREW, A. G. Development and properties of spontaneous oscillations of the membrane potential in inferior olivary neurons in the rat. Dev. Brain Res. 65: 43-50, 1992.[Medline]
• BLEASEL, A. F., PETTIGREW, A. G. The effect of bicarbonate-free artificial cerebrospinal fluid on spontaneous oscillations of the membrane potential in inferior olivary neurons of the rat. Brain Res. 639: 8-20, 1994.[Medline]
• BLOEDEL, J. R., EBNER, T. J. Rhythmic discharge of climbing fibre afferents in response to natural peripheral stimuli in the cat. J. Physiol. Lond. 352: 129-146, 1984.[Abstract/Free Full Text]
• CHRIST, G. J., SPRAY, D. C., BRINK, P. R. Characterization of K currents in cultured human corporal smooth muscle cells. J. Androl. 14: 319-328, 1993.[Abstract/Free Full Text]
• CHRISTIE, M. J., WILLIAMS, J. T., NORTH, R. A. Electrical coupling synchronizes subthreshold activity in locus coeruleus neurons in vitro from neonatal rats. J. Neurosci. 9: 3584-3589, 1989.[Abstract]
• COULTER, D. A., HUGUENARD, J. R., PRINCE, D. A. Calcium currents in rat thalamocortical relay neurones: kinetic properties of the transient, low-threshold current. J. Physiol. Lond. 414: 587-604, 1989.[Abstract/Free Full Text]
• DE ZEEUW, C. I., RUIGROK, T.J.H., HOLSTEGE, J. C., JANSEN, H. G., VOOGD, J. Intracellular labeling of neurons in the medial accessory olive of the cat. II. Ultrastructure of dendritic spines and their GABAergic innervation. J. Comp. Neurol. 300: 478-494, 1990.[Medline]
• DOEDEL, E. J. AUTO: a program for the automatic bifurcation and analysis of autonomous systems. Cong. Num. 30: 265-284, 1981.
• EDELSTEIN-KESHET, L. In: Mathematical Models in Biology., . New York: Random House, 1988
• FITZHUGH, R. Mathematical models for excitation and propagation in nerve. In: Biological Engineering, edited by and H. P. Schwan . New York: McGraw-Hill, 1969, p. 1-85
• GRILLNER, S., MATSUSHIMA, T. The neural network underlying locomotion in lampreysynaptic and cellular mechanisms. Neuron 7: 1-15, 1991.[Medline]
• GUCKENHEIMER, J., LABOURIAU, I. S. Bifurcation of the Hodgkin and Huxley equations: a new twist. Bull. Math. Biol. 55: 937-952, 1993.
• HÁJOS, M., GREENFIELD, S. A. Topographic heterogeneity of substantia nigra neurons: diversity in intrinsic membrane properties and synaptic inputs. Neuroscience 55: 919-934, 1993.[Medline]
• HINES, M. NEURONa program for simulation of nerve equations. In: Neural Systems: Analysis and Modeling, edited by and F. Eeckman . Boston: Kluwer, 1993, p. 127-136
• HOLMES, G. The cerebellum of man. Brain 62: 1-30, 1939.[Free Full Text]
• HOOPER, S. L., MARDER, E. Modulation of the lobster pyloric rhythm by the peptide proctolin. J. Neurosci. 7: 2097-2112, 1987.[Abstract]
• HUGUENARD, J. R., PRINCE, D. A. A novel T-type current underlies prolonged Ca²⁺-dependent burst firing in GABAergic neurons of rat thalamic reticular nucleus. J. Neurosci. 12: 3804-3817, 1992.[Abstract]
• HUTCHEON, B., MIURA, R. M., YAROM, Y., PUIL, E. Low-threshold calcium current and resonance in thalamic neurons: a model of frequency preference. J. Neurophysiol. 71: 583-594, 1994.[Abstract/Free Full Text]
• JACK, J.J.B., NOBLE, D., AND TSIEN, R. W. Electric Current Flow in Excitable Cells. Oxford, UK: Clarendon Press.
• KEATING, J. G., THACH, W. T. Nonclock behavior of inferior olivary neurons: interspike intervals of Purkinje cell complex spike discharge in the awake behaving monkey is random. J. Neurophysiol. 73: 1329-1340, 1995.[Abstract/Free Full Text]
• KEPLER, T. B., MARDER, E., ABBOTT, L. F. The effect of electrical coupling on the frequency of model neuronal oscillators. Science Wash. DC 248: 83-85, 1990.[Abstract/Free Full Text]
• LAMPL, I., YAROM, Y. Subthreshold oscillations of the membrane potential: a functional synchronizing and timing device. J. Neurophysiol. 70: 2181-2186, 1993.[Abstract/Free Full Text]
• LAMPL, I. AND YAROM, Y. Subthreshold oscillations and resonant behavior: two manifestations of the same mechanism. Neuroscience. In press.
• LLINÁS, R. Rebound excitation as the physiological basis for tremor: a biophysical study of the oscillatory properties of mammalian central neurones in vitro. In: Movements Disorders: Tremor, edited by L. J. Findley, and R. Capildeo . London: Macmillan, 1984, p. 135-169
• LLINÁS, R., BAKER, R., SOTELO, C. Electrotonic coupling between neurons in cat inferior olive. J. Neurophysiol. 37: 560-571, 1974.[Free Full Text]
• LLINÁS, R., SASAKI, K. The functional organization of the olivo-cerebellar system as examined by multiple Purkinje cell recordings. Eur. J. Neurosci. 1: 587-602, 1989.[Medline]
• LLINÁS, R., WELSH, J. P. On the cerebellum and motor learning. Curr. Opin. Neurobiol. 3: 958-965, 1993.[Medline]
• LLINÁS, R., YAROM, Y. Electrophysiology of mammalian inferior olivary neurones in vitro. Different types of voltage-dependent ionic conductances. J. Physiol. Lond. 315: 549-567, 1981.[Abstract/Free Full Text]
• LLINÁS, R., YAROM, Y. Oscillatory properties of guinea-pig inferior olivary neurones and their pharmacological modulation: an in vitro study. J. Physiol. Lond. 376: 163-182, 1986a.[Abstract/Free Full Text]
• LLINÁS, R., YAROM, Y. Specific blockage of low threshold calcium channels by high molecular weight alcohols. Soc. Neurosci. Abstr. 12: 174, 1986b.
• MANOR, Y. Construction of an Experimentally-Based Network Model to Explain the Sub-Threshold Oscillations in the Inferior Olive (PhD thesis). Jerusalem: Hebrew University, 1995.
• MARDER, E., ABBOTT, L. F., KEPLER, T. B., HOOPER, S. L. Modification of oscillator function by electrical coupling to nonoscillatory neurons. In: Induced Rhythms in the Brain, edited by E. Basar, and T. H. Bullock . Cambridge, MA: Birkhäuser, 1992, p. 287-296
• MARR, D. A theory of cerebellar cortex. J. Physiol. Lond. 202: 437-470, 1969.[Abstract/Free Full Text]
• MEDA, P., ATWATER, I., GONCALVES, A., BANGHAM, A., ORCI, L., ROJAS, E. The topography of electrical synchrony among -cells in the mouse islet of Langerhans. Q. J. Exp. Physiol. 69: 719-735, 1984.[Medline]
• MEUNIER, C. The electrical coupling of two simple oscillators: load and acceleration effects. Biol. Cybern. 67: 155-164, 1992.[Medline]
• NAYFEH, A. H. In: Perturbation Methods., . New York: Wiley, 1973
• RINZEL, J., ERMENTROUT, B. Analysis of neural excitability and oscillations. In: Methods in Neuronal Modeling. From Synapses to Networks, edited by C. Koch, and I. Segev . Cambridge, MA: MIT Press, 1989, p. 135-169
• ROBINSON, D. A. Adaptive gain control of vestibulo-ocular reflex by the cerebellum. J. Neurophysiol. 39: 954-969, 1976.[Abstract/Free Full Text]
• SASAKI, K., BOWER, J. M., LLINÁS, R. Multiple Purkinje cell recordings in rodent cerebellar cortex. Eur. J. Neurosci. 1: 572-586, 1989.[Medline]
• SHERMAN, A., RINZEL, J. Rhythmogenic effects of weak electrotonic coupling in neuronal models. Proc. Natl. Acad. Sci. USA 89: 2471-2474, 1992.[Abstract/Free Full Text]
• SMOLEN, P., RINZEL, J., SHERMAN, A. Why pancreatic islets burst but single  cells do not. The heterogeneity hypothesis. Biophys. J. 64: 1668-1680, 1993.[Medline]
• SOTELO, C., LLINÁS, R., BAKER, R. Structural study of inferior olivary nucleus of the cat: morphological correlates of electrotonic coupling. J. Neurophysiol. 37: 541-559, 1974.[Free Full Text]
• STROGATZ, S. H. In: Nonlinear Dynamics and Chaos With Applications to Physics, Biology, Chemistry, and Engineering., . Reading, MA: Addison-Wesley, 1994
• WANG, X.-J., RINZEL, J. Alternating and synchronous rhythms in reciprocally inhibitory model neurons. Neural Comput. 4: 84-97, 1992.
• WELSH, J. P., LANG, E. J., SUGIHARA, I., LLINÁS, R. Dynamic organization of motor control within the olivocerebellar system. Nature Lond. 374: 453-457, 1995.[Medline]
• YAROM, Y. Oscillatory behavior of olivary neurons. In: Experimental Brain Research Series. The Olivocerebellar System in Motor Control, edited by and P. Strata . Berlin: Springer-Verlag, 1989, vol. 17, p. 209-220
• YAROM, Y. Rhythmogenesis in a hybrid systeminterconnecting an olivary neuron to an analog network of coupled oscillators. Neuroscience 44: 263-275, 1991.[Medline]

0022-3077/97 $5.00 Copyright ©1997 The American Physiological Society