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