Abstract
Manor, Yair, John Rinzel, Idan Segev, and Yosef Yarom. Lowamplitude 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 singlecompartment, twovariable, HodgkinHuxleylike model for inferior olive neurons. The model consists of a leakage current and a lowthreshold 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 twoparameter 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
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 (sinusoidallike) 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 lowthreshold 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 lowthreshold 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 bicarbonatebuffered 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 linearresistive 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 lowthreshold, 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 lowamplitude oscillations when they are electrically coupled.
METHODS
Extraction of voltageclamp data
The lowthreshold 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 voltageclamp experiments. A detailed description of the experimental results, the spaceclamp problems, experimental protocols, and theoretical validations is given elsewhere (Manor 1995). With voltageclamp data from 15 cells, we determined an activation curve (m
^{3}
_{∞}) and an inactivation curve (h
_{∞}) as shown in Fig. 1
A. The solid curves are the equations
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).
SINGLENEURON MODEL.
Two differential equations govern the dynamics of the single cell
Tools for analyzing the model's dynamic behavior
We are interested in studying the oscillatory (periodic) and steadystate (timeindependent) 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 steadystate 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 twovariable cell model the stability is determined by the eigenvalues of the twodimensional 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, timedependent 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
The steady states and their stability for our system were also explored graphically by phase plane methods (see EdelsteinKeshet 1988 and Strogatz 1994 for general background). A phase plane portrait is a twodimensional plot describing the dynamic relationships between the two dependent variables, in our case V and h;solution trajectories are curves in the Vh plane. The nullclines of the system are the two special curves (not trajectories) along which
In other words, a steady state (for such a twovariable 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 IV 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 (steadystate) 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 twoparameter continuation problem). We used AUTO to explore the behavior of an isolated cell (a twovariable model), or a pair of coupled cells (a fourvariable 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 voltageclamp 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 currentclamp experiments and to construct momentary voltagecurrent (VI) 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 VI curve corresponding to time t′ was constructed by repeating this for several clamping voltages. We note that momentary VI 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
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 steadystate IV 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.
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 ). 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 stimulusresponse 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 (timeindependent) or periodic states, both stable (solid) and unstable (dotted). The curve corresponding to the steady state thus constitutes a cell's steadystate IV relation. A stable cell (A) has a unique stable steadystate V that increases monotonically with applied current. Its IV relation is nearly linear except near the resting potential (I _{app} = 0), where it bends because of the existence of a steadystate (“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 bubbleshaped 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 IV relation develops an Nshaped character, and thus multiple stable steady states may coexist for some current levels. Figure 3 D 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.
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
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 4 A 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. 4 B, the pair oscillates indefinitely for coupling conductances of 0.1 and 0.25 mS/cm^{2}, but not for a coupling conductance of ≥0.5 ms/cm^{2}. 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).
Note that in Fig. 4 C the average channel densities falls in the CO zone (just below the border of the SO zone). Note also that in Fig. 4 D 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 4 E 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^{2}, yet they generate a periodic sustained behavior with a coupling conductance of 0.5 mS/cm^{2}. 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^{2}. A pair was declared as oscillatory if one or two HB points were detected by AUTO. We found no combinations of nonSOs 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 gapjunctionmediated 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 IV 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 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 lowestorder term for the membrane potentials (corresponding to the case of infinite coupling) is V
_{av}, the membrane potential of the average cell. Thus
The coupling current that flows through the gap junctions from cell 2 to cell 1 is defined as
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
Some of these points are illustrated in Figs. 5 and 6 for a case similar to that in Fig. 4 A. We start by describing the case of infinite coupling (Fig. 5 A). 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 Vh 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. 5 B, where the two cells are electrically coupled with g _{coup} = 0.5 mS/cm^{2}.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.
Figure 5 C 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^{2}.
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 _{1}/dt < 0; whenI _{ion,1} < I _{coup}, dV _{1}/dt > 0.
Figure 6 illustrates the accuracy of the strong coupling approximation for I _{coup} derived in Eq. 21. In Fig. 6 A, 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^{2}), the period and minimum amplitude are quite accurate; the maximum amplitude differs by only ∼20%. In Fig. 6 B 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. 6 B we see that when g _{coup} < 0.13 mS/cm^{2}, 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^{2}, the system's fixed point (not shown) is unstable. Figure 6 B 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 conductances—not 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
_{1}/dt is governed mainly by the difference between V
_{1} and V
_{L}, and V
_{1} relaxes quickly to make this difference small, i.e., to achieve V
_{1} ≈ V
_{L}. The other cell (cell 2) is a CO or a conditional bistable cell. The differential equation for its voltage reduces to
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
_{1} ≈ V
_{Ca}. Thus for cell 2 we have approximately
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 VI curves for this case (see Methods). In Fig. 7, momentary VI curves at 3, 10, and 300 ms are plotted for a stable cell (cell 1, Fig. 7 A), a CO (cell 2, Fig. 7 B), an average cell (coupling of cell 1 and cell 2 with an infinite coupling, Fig. 7 C), and cell 1 coupled to cell 2 with a finite coupling (Fig. 7 D). 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 VI 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. 7 C), or the coupled pair of cells (Fig. 7 D), there exists a period of time in which the momentary VI curve has a negative slope at the fixed point (Fig. 7 C, 3 and 10 ms; Fig. 7 D, 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.
DISCUSSION
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 highthreshold 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 Ttype 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 shorttime transient slope conductance is negative. Such a condition may exist only if the voltagedependent 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 lowthreshold calcium current from IO neurons (Fig. 1 A) 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 steadystate 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 twoneuron 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 twocell 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 twoneuron network (Fig. 4 E). 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 stimulusresponse properties for current injection. That is, Fig. 2 distinguishes behavior for ranges of three, not just two, parameters—the 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 voltagedependent 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 twoneuron 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.
Acknowledgments
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 N0001419J1350 and by the United States–Israel Binational Science Foundation (Y. Yarom).
Footnotes

Address for reprint requests: Y. Manor, Volen Center for Complex Systems, Brandeis University, 415 South St., Waltham, MA 02254.
 Copyright © 1997 the American Physiological Society
APPENDIX A: SEGMENTATION OF THE ḡ_{T}g_{L} 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 wedgeshaped region in Fig. 2) is found. With g _{L} fixed, with the use of AUTO we compute the steadystate and periodic solutions as functions of the parameter ḡ _{T}; the results are presented in the bifurcation diagram of Fig. 8 A. We start with ḡ _{T} = 0, where the steady state is V̄= V _{L}, h̄ = h _{∞}(V _{L}). AUTO computes the steadystate 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_{1} and HB_{2}. 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 twoparameter continuation feature. This allows us to track the two HB points as both ḡ _{T} and g _{L} are allowed to change (Fig. 2 B). 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.
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 Nshaped steadystate IV relation, and perhaps be bistable, having two stable steady states for a range of currents. With AUTO, we compute the oneparameter bifurcation diagram: steadystate and periodic solutions as a function of I _{app} (Fig. 9 A). 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_{1} and LP_{2} (where the steadystate solution plot has vertical slope). Now, with twoparameter continuation in AUTO, we trace these points as both I _{app} and ḡ _{T} are allowed to change (Fig. 9 B). 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_{1} and LP_{1} occur are lower than those at which HB_{2} and LP_{2} occur, respectively.
We define four critical values of ḡ _{T}: the value where the two HB points coalesce g _{0}; the two values g _{1}, g _{2} where the HB curve crosses the axis I _{app} = 0; and where the HB_{2} and HB_{1} curves coalesce, g _{3}. Note, the values g _{1}, g _{2} 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 _{0}, at any ḡ _{T} smaller than g _{0} there are no bifurcation points at any injected current. Thus the steadystate solution is unique and always stable, for any I _{app}. Thus the range ḡ _{T} < g _{0} defines the stable neurons.
2) When g _{1} ≤ ḡ _{T} ≤ g _{2}, 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 _{1}, which corresponds to a sliver in the ḡ _{T} g _{L} plane between the ST zone and the SO wedge.
3) When g _{2} < ḡ _{T} < g _{3} (or, in case the HB loop does not cross the ḡ _{T} axis, when g _{0} < ḡ _{T} < g _{3}), bistability cannot occur because I(V _{HB2}) > I(V _{HB1}). g _{2} 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_{1} and HB_{2} at g _{3} marks the lowest ḡ _{T} value for which I(V _{HB2}) ≤ I(V _{HB1}). Thus, for ḡ _{T} ≥ g _{3}, there is a range of negative currents at which two stable steadystate solutions exist. These are the conditional bistable cells.
These previous steps are repeated for different values of g _{L}. The g _{0} values found for different g _{L}s define the border of the ST regime. Note that when g _{0} 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 _{1} and g _{2} as functions of g _{L} define the left and right boundaries, respectively, of the SO wedge in Fig. 2. The g _{3} values found for different g _{L}s mark the border segregating bistable cells from other cells in the ḡ _{T} g _{L} space. Between g _{0} and g _{3}, for a given g _{L}, are located the oscillators, spontaneous and conditional. Their distinction arises after we plot g _{1} and g _{2}.
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. 9 A in which the unstable periodic branch from HB_{2} 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
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 LindstedtPoincaré method, these resonant behaviors are controlled by allowing the unknown frequency of the oscillations, ω_{1}, to also depend on ε (Nayfeh 1973). Thus this frequency is also expanded in series form
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 firstorder terms (σ_{1}, σ̃_{1}, and δ̃_{1}) 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.