|Home | About | Journals | Submit | Contact Us | Français|
Neurons in diverse brain areas can respond to the interruption of a regular stimulus pattern by firing bursts of spikes. Here we describe a simple model which permits such responses to periodic stimuli over a substantial frequency range. Focusing on the omitted stimulus response (OSR) in isolated retinas subjected to periodic patterns of dark flashes, we develop a pharmacologically-based model which accounts for resonances in ON bipolar cells. The bipolar cell terminal contains an LRC oscillator whose inductance is modulated by a transient calcium concentration, thus adjusting its resonant frequency to approximately match that of the stimulus. The model reproduces ganglion cell outputs, which sum the ON and OFF bipolar pathways, and it responds to omitted flashes with approximately constant latencies, as observed experimentally.
Many neural circuits generate responses that embody simple predictions about future events, including cortical networks (Sutton et al., 1967; Klinke et al., 1968; Rogers et al., 1992; Bullock and Karamürsel, 1994; Friedman et al., 2001; Busse and Woldorff, 2003; Jääskeläinen et al., 2004), cerebellar loops (Bullock et al., 1993; Bell et al., 1997), midbrain dopamine systems (Schultz et al., 1997; Schultz, 1998), and subcortical visual pathways (Bullock et al., 1990; Bullock and Karamürsel, 1994; Prechtl and Bullock, 1994; Berry II et al., 1999; Hosoya et al., 2005; Schwartz et al., 2007a; Schwartz et al., 2007b). A particularly simple and vivid example of such pattern recognition occurs when a repetitive sequence is interrupted or ends. Such pattern violations generate characteristic voltage signatures on the human scalp called the omitted stimulus potential or the mismatch negativity (Bullock and Karamürsel, 1994; Klinke et al., 1968), which has been associated with the detection of novelty (Friedman et al., 2001; Rogers et al., 1992; Jääskeläinen et al., 2004). However, the phenomenon does not require higher brain areas or extensive networks: it has been seen in the retina (Bullock et al., 1990; Karamürsel and Bullock, 1994), and recent in vitro studies have shown that it can be robustly and precisely generated by ganglion cells, where it has been called the omitted stimulus response (OSR) (Schwartz et al., 2007a; Schwartz and Berry II, 2008).
Here we focus on the ganglion cell OSR, which has been observed in response to sequences of periodic flashes over the frequency range 6 – 20 Hz, and we investigate a generic mechanism which may be responsible for it. Using a linear model, we show that adaptive conductance dynamics, driven by transient calcium release in ON bipolar cell terminals, can allow a simple peripheral sensory circuit to tune a resonator to respond to interruptions in periodic stimuli over a substantial frequency range.
The paper is organized as follows. We review the OSR in §2, and in §3 describe the resonator mechanism, which builds on standard linear-nonlinear models for the early visual system. An approximate analysis reveals its tuning properties and guides parameter selection, and we confirm that the model reproduces basic experimental phenomena. §4 contains results of simulations and comparisons to experimental OSR data and we summarise and discuss further implications of the work in §5. While certain details are specific to the OSR, we believe that the underlying principles will be generally useful in the study of subthreshold resonant and oscillatory response to periodic stimuli, and to inputs from other oscillatory circuits. Resonance, adaptive tuning and signal cancelation may also contribute to pattern violation detection in other subcortical and cortical areas.
There is a substantial literature on intrinsic oscillatory and resonant properties in both neural circuits and in single neurons. For background on specific ion channels and basic physiological mechanisms, see (Llinas, 1988; Hutcheon and Yarom, 2000; Llinas and Steriade, 2006). For calcium-gated potassium channels similar to those modeled here, see (Llinas and Yarom, 1986). A more detailed Hodgkin-Huxley type model of thalamic neurons, including multiple calcium and calcium-activated currents, is developed in (Hutcheon et al., 1994).
In experiments reported in (Schwartz et al., 2007a; Schwartz and Berry II, 2008), periodic sequences of dark flashes were presented to salamander and mouse retinas in vitro and ganglion cell responses recorded. Sequences were interrupted either by omitting a flash or by terminating the sequence; the OSR occurred in both cases for a subset of the ganglion cells. The population response latency measured with respect to the time of the last flash shifted to later times as the stimulus period increased, but latency relative to the time of the missing flash was constant, suggesting that ganglion cells respond to the omission of a flash from a predictable sequence. Typical cells only produced an OSR over a part of the 6 to 20 Hz range probed in the experiments. At higher frequencies, the timing of the OSR for individual cells shifted continuously as the stimulus period changed, keeping the latency from omitted flash to response constant (Schwartz and Berry II, 2008). At lower frequencies, the timing of the OSR for individual cells shifted in a predictive fashion, but did not always change proportionally with the stimulus period.
Fig. 1 shows typical ganglion cell firing patterns in response to a train of dark flashes with ≈ 50% duty cycle. Those on the left are “start-end” cells that respond strongly to the first and the omitted flashes; those on the right exhibit sustained responses that either decay or grow during the body of the sequence. Although the spike trains of each ganglion cell are reproducibly locked to the stimulus, responses exhibit great variability from cell to cell, e.g., the cell at center right produces essentially no OSR, while that at bottom right rings repeatedly, showing evidence of resonant activity.
Pharmacological studies have shown that blocking the inhibitory effects of amacrine cells does not eliminate the OSR, although it leads to an increase in overall firing rates. However, blocking ON bipolar cell responses using the mGluR6 agonist APB eliminates the OSR, which recovers upon APB washout (Schwartz and Berry II, 2008). This strongly suggests that ON bipolar cells are essential to the phenomenon. Some ganglion cells had a ringing OSR, with two or more peaks of firing, as shown in the lower panels of Fig. 1, and many cells that produced an OSR also exhibited period-doubling during the flash sequence (Schwartz and Berry II, 2008) (not shown here), suggesting that an electrical resonance, perhaps in ON bipolar cells, plays a key role. This led to the development of linear-nonlinear (LN) and resonator bank models (see Appendices D, E and F).
In the next section we build on this work to develop a mathematical model which shows that adaptive conductance dynamics, driven by transient calcium release in ON bipolar cell terminals, can allow a simple sensory circuit to tune itself to respond to interruptions in periodic stimuli over a substantial frequency range. We focus on start-end cells, but note that changes in rectification and relative ON/OFF pathway strengths (e.g., via desensitization) can produce other behaviors. We consider only dark flashes whose duty cycles remain constant at 50% as freqencies vary; bright flashes can also elicit an OSR, but individual ganglion cells respond only to dark or to light flash sequences, not to both. Perturbations in both flash amplitude and interflash durations can also cause an OSR, although again individual cells tend to respond only to one of several perturbation types. See (Schwartz et al., 2007a; Schwartz and Berry II, 2008) for further information on these more subtle effects, which are not addressed here.
Following the LN model (Baccus and Meister, 2002; Dayan and Abbott, 2001), we adopt the linear filters of earlier retinal models and additionally propose that oscillations in the relatively large axonal terminals of ON bipolar cells play a key role in generating the OSR. This assumption is supported by the observation of spontaneous oscillations in bipolar cell terminals of goldfish, in both isolated cells (Burrone and Lagnado, 1997) and retinal slices under light adapted conditions (Protti et al., 2000). An interplay between voltage-gated calcium and calcium-gated potassium channels in isolated cells is responsible for these oscillations. For general background on retinal neuronal circuitry, see (Wässle and Riemann, 1978; Wässle and Boycott, 1991; Cook and Becker, 1995) and (Dayan and Abbott, 2001, Chap. 2).
To formulate the model we choose a linearized version of ion channel dynamics. We make this choice because the retina contains a great diversity of bipolar cells (Burkhardt and Fahey, 1998; Wu et al., 2000; Roska et al., 1998; Roska and Werblin, 2001) containing different kinds of ion channels (Pan, 2000; Mao et al., 2002; Ichinose et al., 2005; Zenisek et al., 2001). Since we do not know which bipolar cells provide input to the ganglion cells that produce an OSR, we would have to estimate many unmeasured parameters in detailed ion channel models of the Hodgkin-Huxley type. A linear model allows us to test essential ideas with minimal free parameters, as in (Richardson et al., 2003).
In this section we develop a model that encompasses retinal input to bipolar cells, their resulting soma and terminal voltages, and output in terms of ganglion cell firing rates.
Bipolar cells do not fire action potentials and so we are concerned with subthreshold oscillations in their membrane potentials. Following standard linear-nonlinear models (Baccus and Meister, 2002; Dayan and Abbott, 2001), we obtain the ON bipolar cell soma voltage VON(t) by convolving the incoming light stimulus s(t) with a linear filter DON(t) representing the cell's temporal receptive field, cf. (Dayan and Abbott, 2001):
Three filters used in our simulations are shown in Fig. 4 below. The kernel functions DON(t) were taken to have a standard biphasic shape with zero mean, thus excluding the effects of slowly changing light level over the whole flash sequence (clearly, the response to on- and offset of each short flash depends on the filter shape). For convenience, we took the difference between two alpha functions: a shape that closely resembles kernels for bipolar cells derived from white noise measurements (Baccus and Meister, 2002). The parameter sets the overall size of bipolar cell voltage fluctuations, but since the circuit described below is essentially a linear oscillator, this simply scales the overall voltage levels thoughout and so we may set = 1 without loss of generality.
The same receptive field with sign reversed, DOFF = −DON(t), is used to compute soma voltages in OFF bipolar cells. These cells have AMPA and kainate-type glutamate receptors in their dendrites that cause their activation levels to desensitize to sustained inputs (DeVries and Schwartz, 1999; DeVries, 2000). To model this we reduce the soma voltage amplitude following the first positive peak at τp, by setting VOFF(t) to 0.7 VOFF(t) for t > τp. We further assume that the OFF cell terminal voltages UOFF(t) are equal to their desensitized soma voltages:
We also examined the effects of mild (10%, 0.9 ) desensitization and no desensitization of OFF cells, obtaining generally similar results but poorer matches to data, as further described in the Discussion.
Voltage dependent conductances, primarily composed of potassium channels, make the terminal dynamics of ON bipolar cells more complicated (Mao et al., 2002; Mao et al., 2003). In general, these may be described by LRC circuits with constant capacitance C and passive conductance gl in parallel (Koch, 1984; Detwiler et al., 1980; Mao et al., 2002), with a third parallel branch containing a voltage-dependent inductance Lk and conductance gk in series (Art and Fettiplace, 1987, Eqn. (2)): see Fig. 2A. The inductance Lk is generally assumed to be proportional to the reciprocal of the derivative of the potassium conductance, which is described by a sigmoidal function of soma voltage, cf. (Jones et al., 1999):
Here Vb is the half-activation voltage at which gk = /2, d is the maximum slope of the sigmoid and the basic timescale is set by and .
We further assume that the LRC circuit receives input current I(t) = βVON(t) proportional to soma voltage, so that its dynamics are described by a parametrically- and additively-forced linear oscillator with time-varying coefficients Lk(t), gk(t) that depend on soma voltage:
Here UON(t) describes the ON cell terminal voltage which, combined with the output UOFF(t) of OFF bipolar cells, determines ganglion cell input via glutamate release. Details of the derivation of (4) are provided in Appendix A, and a circuit diagram and illustration of the impulse response function appears in Fig. 2 (panels A and B).
Terminal voltages of both ON and OFF bipolar cells govern glutamate release and effectively determine the firing rates of the postsynaptic ganglion cells. For simplicity, we assume a linear relationship between terminal voltage and glutamate release in the range of interest, and a rectified linear relationship between glutamate release and ganglion cell firing rates f(t), so that the net effect is described by
where U(t) = UON(t) + UOFF(t) denotes the (summed) terminal voltage and Uθ is the rectification threshold. Like in Eqn. (1), is an output scale factor, and adjustment of it and of the threshold Uθ allows one to match firing rates of specific ganglion cells and to encompass both start-end and sustained firing cell types. While this is probably adequate for ON bipolar cells, which are observed to operate around a voltage above the rectification threshold, rectification of UOFF(t) prior to summation may be more realistic (Demb et al., 2001; Margolis and Detwiler, 2007). We tested this situation in §4.
Due to the coupled dynamics of voltage-gated calcium and calcium-gated potassium channels, calcium concentration substantially influences the potassium conductances in bipolar cell terminals. Specifically, increases in calcium concentration shift the sigmoidal potassium conductance vs voltage curve of Eqns. (3) to the left, cf. (Jones et al., 1999, Fig. 4). To incorporate this effect we set Vb = Vb0 − α, where is the calcium concentration determined by VON, Vb0 is the half-activation voltage at zero calcium concentration and α is a slope factor which also sets the units of the last term. Logarithmic dependence on was also tried in place of this linear function, without significant differences in model predictions (results not shown).
Detailed studies of calcium-gated potassium conductances have not been carried out for salamander bipolar cells, so we follow the work of Art, Fettiplace et al. (Art and Fettiplace, 1987; Jones et al., 1999) and assume that the effect of calcium on the potassium conductance dominates that of soma voltage variations. Specifically, writing , Eqns. (3) become
where we take the offset b = (Vb0 − VON)/α as constant and α is absorbed into the slope d. Together, b and d describe how calcium tunes conductance.
Since calcium channels are activated by depolarization, calcium concentration (t) can be described by a leaky integrator driven by the rectified soma voltage:
where βv can also be set equal to 1 since it simply sets the calcium level, which is already scaled by d in Eqns. (6). With appropriate time constants, τCa, the resulting concentrations fluctuate mildly as they rise and fall, while their time-averaged levels over the stimulus duration differ for different stimulus frequencies, thereby tuning the resonant frequency of the terminal. See Fig. 3. The removal of VON-dependence in Lk and gk and the relatively slow calcium dynamics (t) imply that Eqn. (4) behaves like a classical forced and damped linear oscillator: a fact that we appeal to in the analysis of §3.3 and parameter discussion of §3.4.
Fig. 2 summarizes the model, which adds a calcium-tuned terminal oscillator to the ON bipolar cell pathway of the LN model. In the ON branch: 1) the stimulus enters through a bandpass filter to produce the soma voltage; 2) the rectified soma voltage then goes through a low pass filter to produce calcium concentration in the ON cell terminal; 3) calcium concentration adjusts the variable inductance and resistance branch of the LRC circuit as the soma voltage waxes and wanes, and 4) current input proportional to soma voltage passes through the LRC circuit to produce the output terminal voltage. 5) The OFF branch uses the ON bandpass filter but with sign reversed and additional desensitization, and it lacks the terminal dynamics. 6) Finally, the summed outputs from ON and OFF pathways are rectified to produce the ganglion cell firing rate.
If calcium concentration were constant, the terminal conductance gk and inductance Lk would also be constant, along with C and gl, and in the case of subcritical damping, the resonant frequency of the ON terminal circuit would be given by:
(using Eqns. (6)). A brief calculation given in Appendix B reveals that the frequency vs. calcium concentration curve is monotonically decreasing for > b. Thus, stimulus sequences with lower frequencies that release more calcium produce lower resonant frequencies in the ON terminals. Hence calcium can provide a self-tuning mechanism to match the terminal resonance with the stimulus, and, as emphasized in §3.4, once is determined, just four parameter groups describe this process: C, gl/, d and b. The first two are analogous to the conductance-capacitance ratios that determine resonant frequencies and bandwidths in (Richardson et al., 2003).
As we have noted, calcium concentrations (t) fluctuate along with terminal voltages during stimulus presentation, and lower stimulus frequencies cause larger fluctuations, due to higher and wider rectified soma voltage peaks. However, for sufficiently long sequences (~ 10 or more flashes), the instantaneous frequency f0(t), computed from Eqn. (8) using (t), remains approximately constant during the latter part of the sequence, as illustrated in Fig. 3. To estimate the terminal frequency we therefore define a shortterm average as follows:
where T− < Tmax < T+ are defined by setting (T+) = (T−) = 0.75 (Tmax), (t) reaches its maximum value at t = Tmax, and T− and T+ are the first times during the stimulus presentation at which the 3/4-maximum value is reached as (t) rises and falls (middle panel of Fig. 3).
For the receptive field of Fig. 2 and calcium time constant τCa = 300 ms, we find that varies with stimulus frequency as shown in the upper middle panel of Fig. 4, being approximately parabolic with a peak at 11 Hz. Our analysis therefore implies that falling values of beyond the peak will cause the resonant frequency (8) to rise, suggesting that the circuit may be able to accomodate to stimulus frequencies over a range above 11 Hz. After some experimentation, we found parameters that provide a reasonable match over the range 12 – 18 Hz, as shown in the upper right panel of Fig. 4 ( = 15 Hz/mV, Uθ = 35 mV, C = 2 pF, = 4.3 MH, d = 0.1, b = 9.5, k = gl/ = 0.0025 (gl = 0.01 nS, = 4 nS) and β = 28pA/V). We refer to these as the high frequency parameter set.
Since we appeal to the falling branch of the vs frequency curve to the right of the peak, slower receptive field dynamics are required to cover the lower stimulus frequency range. Taking the filter shown in the lower left panel of Fig. 4, with all other parameters as above except d = 0.06 and b = 13 (the low frequency parameter set), the peak in calcium concentration occurs near 8 Hz and we obtain the match over the range 8 – 13 Hz shown in black in the lower right panel of Fig. 4. The fact that this is achieved by adjusting b and d predicts that the properties and density of calcium channels in cell terminals should predominate in determining ON cell contributions to OSR. A slightly slower filter produces a better frequency match (green curves), but provides poorer fits to latency data for ganglion cells (see Fig 7).
The quasisteady approximation for calcium concentration is used to guide parameter selection only in this section; all the simulations that follow were obtained by integrating the full model equations.
In addition to the filter kernel DON(t) and desensitization in the OFF channel, set at 0.7, full specification of the model requires 12 parameters:
but only three are critical in obtaining the OSR. The scale factors , βv and can all be set to unity without changing the results. (As noted above, βv is absorbed into d, and and merely set bipolar cell soma voltage and ganglion cell firing levels, all of which are stages in an essentially linear system.) Desensitization, the LRC input current factor β, and the ganglion cell threshold Uθ do not affect OSR latency, although they change the form of the response, allowing us to reproduce both start-end and sustained firing cells. The calcium release time constant τCa is set at 300 ms, but values throughout the range 275 – 500 ms yield similar results.
The remaining six parameters determine the natural frequency of the LRC oscillator, but Eqn. (8) shows that it depends on four quantities: the timescale , the slope d and offset b of the sigmoid describing how calcium tunes the frequency, and the conductance ratio k = gl/. Numerical experiments show that, provided k 1 (consistent with subcritical damping) and the filter timescale is appropriate, only the first three of these require careful tuning to obtain realistic OSRs (indeed, only d and b are changed in passing from the high to the low frequency parameter set). In summary, only three of the twelve parameters listed in (10) are critical to the final results, and the model is more robust than its multiparameter form might suggest. Further details can be obtained from the first author.
To further test the model, we computed spike-triggered averages (STAs) from the ganglion cell firing rate response f(t) to white noise inputs (Dayan and Abbott, 2001). Fig. 5 (left panels) confirms that STAs obtained with both parameter sets are consistent with typical receptive field dynamics of ganglion cells observed in white noise experiments. The ON pathway bipolar terminal dynamics adds an oscillatory tail to the biphasic response of the filter alone, making the STA weakly triphasic, but this shape is consistent with STAs for many ganglion cells measured at the higher light levels used in previous studies of the OSR. See Appendix C for examples, and compare with biphasic STAs for lower light intensities in (Baccus and Meister, 2002; Segev et al., 2006).
We also examined the model's response to step changes in intensity, confirming that a short flash evokes one burst of firing, and two bursts follow the onset and offset of a longer dark step, as shown in Fig. 6. These firing patterns resemble those of most ON/OFF-type ganglion cells, which comprise up to 75% of all ganglion cells in the salamander (Burkhardt et al., 1998; Segev et al., 2006; Schwartz et al., 2007a). (The ON pathway output had to be increased by adjusting β to bring it closer to that of the OFF pathway to obtain the second burst for the long flash; this is reasonable since the response to a long flash increases with increasing flash duration, indicating that light adaptation is involved in increasing the gain of the ON response.)
To simulate the effect of the drug APB, which blocks light responses in ON bipolar cells, we also computed the response of the model with the ON pathway blocked (UON(t) 0). As shown in Fig. 5 (right panels), there is a substantial response to the first flash, little activity thereafter, and no OSR, closely resembling the experimental data (Schwartz and Berry II, 2008, Fig. 8b).
To confirm that the calcium dynamics and terminal resonator mechanism is required for the two-pathway model to produce a realistic OSR, we carried out simulations of a simple LN model, using the fast filter of Fig. 4, but with the terminal dynamics disabled. Here the input to the ganglion cells is simply the sum VON + VOFF. We found that no OSR is produced in this case: an N-flash sequence elicits N peaks in the ganglion cell response, the final one being the “direct” response to the last flash. Indeed, while this final peak can occur after the omitted flash for higher frequency stimuli, the latency between the start of the omitted flash and this peak varies from −5 to +10 ms over the frequency range 12 – 18 Hz, i.e., at 12 Hz the final peak precedes the omitted flash. These results are not changed by choices of parameters such as the desensitization. The experimental data of Fig 7 show that latencies lie in the range of 66 – 84 ms for 12 – 16Hz, and remain nearly constant for a given cell. A sample of LN model data in shown in Appendix D. We also investigated two second order LN models as described in Appendix F and found that they, too, cannot produce the nearly-constant latency from omitted flash to OSR.
We now describe the model's response to sequences of dark flashes. We remark that duty cycles were kept constant at ≈50% as the frequency varied in some of the experiments reported in (Schwartz et al., 2007a; Schwartz and Berry II, 2008), while others were done with constant flash duration. For the model we use a 50% duty cycle throughout. Samples of results for the high and low frequency parameter sets appear in Fig. 7 (left and right columns, respectively). The bipolar cell outputs UON and UOFF and the summed ganglion cell voltage trace with firing threshold Uθ illustrate how the components of the model of Fig. 2 combine to produce the OSR.
The analysis of §3.3 suggests that specific bipolar cell parameters can produce constant latency OSRs over only part of the relevant frequency range. The simulations confirm this. With the fast receptive field filter of Fig. 4 and the high frequency parameter set, robust OSRs are obtained for stimuli in the range 12 – 18 Hz. The bottom left panel of Fig. 7 shows that latencies from the start of the omitted flash (vertical green bars) to the peak in firing rate of the OSR remain approximately constant (triangles), varying only from 74 to 83 ms. As in (Schwartz et al., 2007a; Schwartz and Berry II, 2008) we plot latencies vs. stimulus period. Very similar behavior was seen for three ganglion cells that were tested over the same range (dashed black lines; latencies measured with respect to the time of the last flash). The data indicate considerable variability among cells (≈ 10 ms) in absolute latency, but the slopes of last flash to OSR curves are reproduced accurately by the model.
The right column of Fig. 7 shows model simulations using the low frequency parameter set, also subject to a 12 Hz stimulus, as well as over the range 8 – 13 Hz (bottom). Again a robust OSR is observed although it is smaller than the response to the first flash. Absolute latencies for this slower model cell are significantly longer than those for the fast cell, since it has slower filters for the bipolar cells. The delay between the last flash and the OSR (circles in bottom panel) shifts in a predictive direction as the stimulus period increases. However, latencies measured with respect to the omitted flash are not constant (triangles). Similar behavior was observed in the experimental data (dashed black) over this frequency range. Model cells with different parameters can produce OSR latencies that are more nearly constant with respect to the omitted flash (cf. green curves in Fig. 4), but here we present model results parameterized to better match the data. The model also captures the experimental observations that both the fraction of cells exhibiting an OSR, and OSR strengths, decline with frequency. As noted above, calcium concentrations fluctuate more widely for low frequency stimuli, implying greater time variability in f0(t) and poorer matching between the circuit's resonant frequency and stimulus frequency.
In comparing model results and data in Fig. 7, one should recall that typical ganglion cells each produce OSRs over only part of the 6 – 20 Hz frequency range that was investigated experimentally. The data of (Schwartz et al., 2007a), including its Supplementary Materials, and of (Schwartz and Berry II, 2008, Figs. 1 and and3)3) (cf. Fig. 1), also indicate considerable variability in both OSR intensity and latencies from cell to cell, as reflected in the range of absolute latencies in the experimental data of Fig. 7. The model can reproduce this by using filters, DON(t) and DOFF(t), with different times to peak, but for clarity we present only a single choice. Similarly, peak model firing rates are between 50 and 150 Hz, similar to the 100 – 200 Hz data of (Schwartz et al., 2007a), but rates in the range 20 – 40 Hz have also been observed, as illustrated in the lower left panels of Fig. 1. The free parameter that scales the ganglion cell firing rate allows the model to encompass this variability; here it was kept fixed for all simulations.
As shown in Fig. 7, the OSR arises from a subtle interaction between the resonant oscillatory response of the ON bipolar cell terminals and the linearly filtered and desensitized OFF bipolar cell response. Our model does not rectify either signal until after they are summed (Eqn. (5)), but this simplification is not necessary. In Fig. 8 we show that an OSR is also produced when the OFF pathway output is rectified prior to summation with UON(t). While the OSR is relatively weak, we note that the rising firing rates prior to the last flash are qualitatively similar to responses of facilitation type cells: see Fig. 1 (bottom right) and the facilitating and weak OSR cell of (Schwartz and Berry II, 2008, Fig. 13B). Thus, differences among the properties of bipolar cells feeding into a given ganglion cell could help explain the diversity of ganglion cell response types observed experimentally.
In this paper we show that a simple LRC oscillator, with inductive and resistive elements modulated by release of a relatively slow ionic species, can dynamically adapt its resonant frequency to that of a periodic stimulus. In combination with other, static, circuit elements, this plays a central role in producing a response that follows an omitted event by a fixed latency. An earlier study of linearized conductance-based models also found that subthreshold resonances can be tuned by conductance and capacitance parameters (Richardson et al., 2003), although that oscillatory mechanism involves interactions between a single voltage equation and gating variables, while ours derives from membrane capacitance and inductance.
Our dual pathway model, including standard OFF bipolar cells and ON bipolar cells whose terminals adaptively resonate at the stimulus frequency, accounts for the omitted stimulus response (OSR) observed in isolated mouse and salamander retinas. The resonance arises from the interaction between voltage-gated calcium channels and calcium-gated potassium channels, as described in bipolar cells (Burrone and Lagnado, 1997; Protti et al., 2000) and other neurons (Art and Fettiplace, 1987). Unlike previous models, the electrical resonance adapts its frequency using calcium-mediated feedback. Such a mechanism is expected to operate for any resonator involving BK-type potassium channels, which are gated both by voltage and calcium (Hille, 1992; Jones et al., 1999). Thus, the key novel feature relies on well-known properties of ion channels in bipolar cells (Burrone and Lagnado, 1997; Sakaba et al., 1997). Moreover, since the mechanism is basically linear, and does not rely on subtle current-voltage characteristics, it may be more generally responsible for predictive responses to disruptions in regular patterns such as those noted in the Introduction.
The model captures several qualitative features of the experimental observations. First, it produces OSRs with constant latencies from the expected flash for the frequency range 12 – 18 Hz (Fig. 7, upper panel). Second, over the lower frequency range, 8 – 13 Hz, latencies of model OSRs shift in the predictive direction, albeit not proportionally to the stimulus period (Fig. 7, lower panel), and OSR strengths are relatively smaller and less robust (robustness results not shown). Third, with fixed parameters the model can produce an OSR either for the high or low frequency range, but not both. All these behaviors match experimental data from individual ganglion cells. Finally, the spike-triggered averages and responses to flash sequences with the ON bipolar cells blocked match those observed for ganglion cells, and the OSR survives rectification in the OFF bipolar pathway prior to summation at the ganglion cell.
As described in Appendix E, we also explored an alternative mechanism in which ganglion cells receive input from a spectrum of resonating bipolar cells with differing, but fixed, resonant frequencies. This resonator bank model generally fails to produce OSRs over a realistic frequency range. The addition of gap junction coupling can improve the range, but without producing constant latencies from expected flashes. This work does not rule out these or similar fixed parameter models, but it does suggest an important role for adaptation. It also shows that changes in stimuli can cause adaptation via calcium influx, with the natural prediction that manipulations that clamp intracellular calcium concentrations would lead to OSRs that cannot maintain a constant latency.
Our model includes desensitization in the OFF bipolar pathway, as has been observed experimentally to arise from the properties of ionotropic glutamate receptors (DeVries, 2000). Because different degrees of desensitization are found for kainate- versus AMPA-type glutamate (DeVries and Schwartz, 1999), we tried different degrees of desensitization in the model and found that the existence and timing of OSR was unaffected, even with no desensitization at all (results not shown). However, neglecting desensitization in the OFF pathway can lead to responses to second and succeeding flashes that are relatively large compared to the first flash response: behavior that is rarely observed. In addition, desensitization often caused the OFF bipolar pathway, alone, to fail to exceed threshold after the first flash (simulations not shown). This property is important to match the effects of the drug APB, which blocks ON bipolar cells, abolishing the OSR and the sustained response (Schwartz and Berry II, 2008).
In addition to bipolar cell terminals, adaptation can occur in photoreceptors, bipolar cell somas and ganglion cells. Our model ignores such detailed dynamics by describing the initial processing from photoreceptors to bipolar cell somas by a simple linear filter, and using a rectified linear relationship between terminal voltage and ganglion cell firing rate to produce the final output. Inclusion of these dynamics was not necessary to explain the salient features of the OSR, although they might be important to correctly predict firing rates of individual ganglion cells. The current model cannot produce OSRs for the more complex pattern violations observed in (Schwartz et al., 2007a; Schwartz and Berry II, 2008); nor does it produce the period doubling and higher harmonic patterns shown there. Addition of adaptive dynamics in other cells might solve this problem.
The adaptive resonator model was motivated specifically by the omitted stimulus response in the retina (Schwartz et al., 2007a), but its feedback mechanism, in which calcium accumulation changes the effective inductance of calcium-gated potassium channels, may be quite generic. This mechanism relies on well-established properties of BK-type ion channels, which are found widely in the brain. Thus, the ability of a resonator to adapt its characteristic frequency over some range of input frequencies might be active in cellular oscillators in other subcortical and cortical areas (Llinas and Yarom, 1986; Llinas, 1988; Puil et al., 1994; Gutfreud et al., 1995; Hutcheon et al., 1996; Pike et al., 2000; Hutcheon et al., 1994). An important factor is the accumulation of calcium with a time constant longer than each cycle of the stimulus, which allows the steady-state calcium level to track the stimulus frequency (Fig. 3). Cells with faster calcium dynamics will not significantly accumulate calcium, allowing resonators to be either fixed or adaptive depending on the calcium pumps and buffering.
Similar feedback may also be possible using the properties of Ih channels (Robinson and Siegelbaum, 2003). These channels are activated by hyperpolarization and produce a depolarizing current, which can drive a resonance. They also typically have long time constants, and so might encode stimulus frequency in a manner analogous to calcium accumulation in the current model. Ih channels exhibit diverse expression patterns in bipolar cells, with mixtures of all four molecular types (HCN1-4) in different bipolar cells (Ivanova and Müller, 2006), as well as in the axon terminals of ON bipolars (Fyk-Kolodziej and Pourcho, 2007). Interestingly, Ih channels can be gated by cyclic nucleotides in addition to voltage (Craven and Zagotta, 2006). Because the metabotropic glutamate receptors on ON bipolar cells act via changes in cGMP levels, Ih channels may interact in a complex way with this signaling cascade. Since our calcium feedback mechanism successfully captures the OSR observations, we have not attempted to add Ih channels to our model.
J.G. and P.H. were partially supported by PHS grants MH58480 and MH62196 (Cognitive and Neural Mechanisms of Conflict and Control, Silvio M. Conte Center). G.S. and M.B. were supported by grants EY14196 and EY17934 from the NEI. We thank the referees for valuable comments.
Denoting the voltage across the circuit of Fig. 2(A) by UON(t), the currents through the first two channels are described by
where C and gl are the capacitance and leak conductances. The current through the third, variable inductance and conductance channel is given by
Assuming that the cell is in its resting state at stimulus onset t = 0, we have I3(0) = 0 and (12) may be solved to yield
where and . Summing the currents in the three channels and using (11) and (13), we find that the dynamics governing the terminal voltage is described by
where I(t) = I1(t) + I2(t) + I3(t) is the total current flowing through the terminal. Multiplying both sides of (14) by and taking derivatives with respect to t yields Eqn. (4).
Here we verify that the frequency vs. calcium concentration curve is monotonically decreasing for > b. Noting that the sigmoid expression 1 + e−4d(−b) in the denominator of Eqn. (8) is a monotonically decreasing function of , we define Φ = 1 + e−4d(−b) and rewrite Eqn. (8) in terms of this new variable:
where we also note that k = gl/ (0, 1) and Φ (1, ∞). Computing the derivative of this expression with respect to Φ:
we find that it is always positive for Φ < 2 or, equivalently, for > b. We conclude that, as increases, both Φ and f0 decrease.
In the experimental studies on which this paper draws, visual stimulation was provided by a light-emitting diode (LED). However, many previous studies that have displayed STAs for retinal ganglion cells in the salamander used a computer monitor to generate visual stimuli (Baccus and Meister, 2002; M. Meister and Baylor, 1994); (Segev et al., 2006, Figs. 1 and and2).2). The experiments of (Schwartz et al., 2007a; Schwartz and Berry II, 2008) differ in that the LED stimulation had a higher average light level than was possible with a computer monitor. Moreover, some previous studies used flickering checkerboard stimuli to generate the STA and displayed the temporal profile of the receptive field center (Baccus and Meister, 2002; M. Meister and Baylor, 1994; Segev et al., 2006). This is not equivalent to the STA under spatially-uniform flicker, because it excludes the influence of the receptive field surround. We therefore display here typical examples of STAs taken with the same LED used in (Schwartz et al., 2007a; Schwartz and Berry II, 2008), for direct comparison with the model filters shown in Fig. 4.
As seen previously (Segev et al., 2006), there was considerable diversity in STAs among different ganglion cells. Unlike previous studies, we found that under LED stimulation many cells had significantly triphasic STAs. Fig. 9 shows multiple examples from three different subtypes of triphasic STAs (blue: strong triphasic; purple: fast weak triphasic; green: slow weak triphasic). As described previously (Segev et al., 2006), these groupings should not be viewed as discrete functional classes, but instead as representatives of a continuum of functional properties. Also shown for comparison are examples of purely biphasic STAs (red). Finally, note that the latency of the first peak in the STAs typically is considerably shorter than in previous studies (30-50 ms), presumably because of the higher light level from the LED.
Cell types were distributed as follows in this data set: triphasic STAs (13/43 cells), biphasic STAs (10/43 cells), monophasic STAs (8/43 cells), other response types (12/43 cells).
In Fig. 10 we show an example of the response of the LN model. As noted in §3.5, no OSR occurs: each dark flash produces one burst of spikes, due to desensitization in the OFF pathway. Moreover, in the 12 Hz case shown on the left panel the final peak precedes the missing flash. Response latencies vary significantly and remain substantially smaller than those of the experimental data. Compare with Fig. 7 (top panels). Further LN models are discussed in Appendix F.
It is well known that in the retina multiple bipolar cells map to the same ganglion cell. In the adaptive resonator model, we assume that these bipolar cells have similar properties, hence the output of all ON bipolar cell terminals is represented by a single terminal. In contrast, the resonator bank model emphasizes differences among the properties of bipolar cell terminals. As in the adaptive resonator model, signal processing from photoreceptors to all ON bipolar cell somas is assumed to be identical, and the terminals are described as LRC circuits, but they are endowed with a range of resonant frequencies. Different phase shifts in the resulting terminal responses cause reinforcement and cancellation in the summed input to the ganglion cell that can produce OSRs.
Specifically, we ignore calcium-dependent ion channel dynamics and assume that each terminal is a damped linear LRC oscillator with common conductances gl, gk and capacitance C, and a unique inductance value , implying a unique resonant frequency. The LRC circuitry for each terminal is the same as in the adaptive resonator model of Fig. 2(A) and the dynamics for the ith terminal is also described by
Eqn. (17) differs from the adaptive resonator model (Eqn. (4) above) in that the inductance and conductance are time-independent, and the inductance differs for every terminal. As in the main text we assume a linear relationship between terminal voltages and glutamate release, and a rectified linear relationship between glutamate release and ganglion cell firing rates. This leads to the same expression for ganglion cell firing rates f(t) as in Eqn. (3) of the main text, but with .
With appropriate parameter values the resonator bank model can produce OSRs, but only over a limited frequency range. Consider a model whose ON bipolar cell terminals have natural frequencies spanning the range 12 – 18 Hz. With the faster receptive field shown in Fig. 4 of the main text and the high frequency parameter set (except Uθ = 10 mV and β = 8pA/V in order to scale the contribution of the seven terminals to have approximately the same weight as before for the entire ON pathway), one can produce OSRs for stimuli in the range 12 – 15 Hz. However for 12 and 13 Hz, spikes persist throughout the flash sequences. Furthermore, no OSR occurs for 16 – 17 Hz, and two very weak OSR peaks are seen at 18 Hz. In all cases, OSR sizes are much smaller than the first firing rate peak, and in some cases, smaller than those in the middle of the sequence (Fig. 11A). Similar results are seen for the low frequency parameter set. This lack of robustness appears intrinsic, as we now explain.
Responses of different terminal oscillators exhibit different phase shifts in relation to the flashes. These phase shifts depend on stimulus period and they accumulate over time, leading to variable peak amplitudes and positions in the summed response (see Fig. 11B). Although the summed response is most strongly influenced by the oscillator tuned to the stimulus frequency, subtle interactions among the phases of other oscillators cause the OSR size and latency to vary as the number of flashes changes (results not shown), even for stimuli in a narrow frequency range. Lowering the damping term produces larger amplitude OSRs and more stable latencies from expected flash to OSR, but it also leads to multiple OSR peaks. A damping rate that is very low when the stimulus is on and high when it is off yields improved results, but at the expense of another form of adaptation.
Anatomical studies reveal widespread gap junctions among bipolar cells (Saito and Kujiraoka, 1988; Umino et al., 1993; Cook and Becker, 1995), which are generally believed to synchronize coupled cells. We tested the effects of gap junctions and synchronization for the parameter set specified in previous section. Gap junctions are modeled as linear ionic conductors connecting different LRC circuits. Assuming all-to-all connections in a network of N ON bipolar cells, the input current on the right hand side of Eqn. (14) is supplemented by the term , so that the dynamics of the ith terminal is now described by
where the coupling strength g is taken to be the same between different cells, and the superscript j denotes the jth ON bipolar cell terminal.
Numerical experiments show that with coupling strengths ≥ 0.05 nS, the voltages of different terminals synchronize, implying fixed phase shifts among the oscillators (Fig. 12A). The phase of the summed responses therefore approximately matches that of the oscillator whose resonant frequency lies in the middle of the range (15 Hz here). With all terminal responses synchronized, more reasonable OSRs can be obtained, albeit not over the entire frequency range (Fig. 12B). Latencies from expected flash to the peak OSR vary widely (from 53 to 85 milliseconds), in contrast to the near-constant latencies from last flashes (Fig. 12C), and this behavior is not restricted to a specific parameter set. It could possibly be improved by introducing a feedback loop which adjusts response phases, but this would again involve adding a further form of adaptation.
There is a substantial literature on models of ganglion cell function. Starting with the early work of Rodieck, quantitative models used the spatial and temporal dynamics of the receptive field to predict the time-varying firing rate of a ganglion cell (Rodieck and Stone, 1965a; Rodieck and Stone, 1965b). Such models involve convolution of the stimulus with a linear filter, followed by rectification. More general nonlinear functions have been used, resulting in linear-nonlinear (LN) models (Chichilnisky, 2001), in which the linear kernel describes processing carried out by the classical receptive field. These models have been extended to second order, using Wiener kernel methods (Marmarelis and Marmarelis, 1978; Sakai and Naka, 1987), covariance analysis (Fairhall et al., 2006; Pillow and Simoncelli, 2006), and a dual-LN model of ON/OFF cells (Geffen et al., 2007). Other models have predicted individual spike times using a linear filter followed by a threshold and spike afterpotential (Keat et al., 2001; Pillow et al., 2005), or have included refractoriness (Berry II and Meister, 1998).
Naturally, one wonders if such models, especially second-order ones, can explain the OSR. To this end, we explored two forms of second-order model: covariance and dual-LN. Covariance models have been shown to explain up to 100% of the information contained in spike times under white noise stimulation (Fairhall et al., 2006). The dual-LN model is often similar in practice to the covariance model and allows the interpretation that its two terms represent processing in the ON and OFF channels, mediated primarily through ON and OFF bipolar cells (Geffen et al., 2007). In each case, we used 30-60 min recordings of ganglion cell spiking under spatially uniform, temporally modulated white noise to derive a model.
The covariance model was formed as described in (Fairhall et al., 2006). In brief, the matrix of light intensity covariations at all pairs of times before a spike was compiled over all occurring spikes. This matrix was diagonalized, and the eigenvalue spectrum examined to determine how many stimulus features (each a corresponding eigenmode) were significant. Typically, two such features were identified. To formulate the dual LN model, the spike-triggered stimulus average (STA) was compiled over all occurring spikes (Geffen et al., 2007). Typically, the total STA was OFF-type. Then for each spike, the stimulus was convolved with the total STA. If the resulting value was positive (for an OFF-type total STA), the spike was classified as “OFF-type” and if it was negative, as “ON-type.” Finally, we calculated separate ON- and OFF-type STAs using the respective subsets of spikes.
We found that neither model could account for the predictive timing of the OSR. In both cases, the peak response locked approximately to the time of the last flash, rather than to the time at which the next flash would have occurred.
Fig. 13 shows an example cell analyzed with the covariance model. This cell has two significant eigenmodes; one of ON-type and the other OFF-type (Fig. 13A). Note that the latter has a shorter latency than the ON-type mode and also that it is somewhat triphasic. Each mode was convolved with a set of periodic stimuli having 12 flashes, 50% duty cycle, and different frequencies, and the response was aligned to the time of the last flash (Fig. 13C). We found that the ON-type mode produced a final peak of excitation, which was roughly invariant with respect to the time of the last ON flash, with a latency of ~115 ms (Fig. 13D). The OFF-type mode produced an initial trough of inhibition, followed by excitation (Fig. 13E), again roughly invariant with respect to the last flash (Fig. 13F). Examining the timing of the peak in excitation from ON and OFF modes, we found that both showed a slight decrease in latency as the stimulus period increased (Fig. 13B, circles). This occurred because longer periods matched the bandwidth of the filters better and hence had a somewhat higher effective contrast. This behavior is strikingly different from a constant delay with respect to the time of the expected flash, as found for the OSR (Fig. 13B, dotted).
Fig. 14 shows analogous results for the dual LN model. Here, the OFF STA again has a shorter latency than the ON STA, but it is not appreciably triphasic (Fig. 14A). When the ON STA was convolved with stimuli at different frequencies, it produced a peak of excitation that was roughly invariant with respect to the last flash (Figs. 14C,D). The OFF STA produced only inhibition at the end of the flash sequence (Figs. 14E,F). The latency of excitation from the ON STA again decreased slightly as the stimulus period increased, in striking disagreement with the latency of the OSR (Fig. 14B).
Both the covariance and dual LN models involve an additional nonlinearity. We have not implemented this step here, because a single static nonlinearity clearly cannot shift the OSR peak by different amounts for different frequencies. Nor have we implemented mechanisms of contrast gain control (Shapley and Victor, 1979; Berry II et al., 1999), because they also would not shift the peak differentially. Results are shown here for two example cells, but we observed no ganglion cells with behavior qualitatively different from these.
More recently, (B. Werner and Passaglia, 2008) have studied the OSR and concluded that the predictive shift in response timing can be explained by the ON channel of a dual LN model. Unfortunately, this result rests on a mistaken definition of response latency. In previous experimental studies (Schwartz et al., 2007a; Schwartz and Berry II, 2008), latency was measured with respect to the peak in the firing rate of a ganglion cell, while in (B. Werner and Passaglia, 2008), it is measured as the onset of the input current to a ganglion cell. This is a crucial difference. In an LN model, the onset of firing can be adjusted by changing the threshold of the static nonlinear function, while the peak time does not change, regardless of the form of nonlinear function. In fact, this is precisely how the “timing” of the OSR changed in the model of (B. Werner and Passaglia, 2008).
In Fig. 15 we demonstrate this property for the same LN models studied previously using a stimulus with flashes of constant 40 ms duration, as in (B. Werner and Passaglia, 2008). If we apply a threshold to the output of the ON channel, then this onset time can increase as the stimulus period increases, for at least part of that range. The exact form of this shift depends sensitively on the details of the linear filter (cf. Figs. 15A and 15B). Furthermore, a roughly proportional increase of latency versus stimulus period can sometimes occur even if the absolute latency is not predictive of the time of the next flash (cf. Figs. 15C and 15D). Note that for the same models, the time of the peak output from the ON channel does not shift in a predictive fashion (Figs. 13 and and14).14). We emphasize again that the experimental data shows that the time of peak firing shifts predictively as a function of the stimulus period, and hence, this is the property that any model must explain. Thus, we conclude that neither of these two types of second-order model can explain the OSR.