## Abstract

A key problem in computational neuroscience is to find simple, tractable models that are nevertheless flexible enough to capture the response properties of real neurons. Here we examine the capabilities of recurrent point process models known as Poisson generalized linear models (GLMs). These models are defined by a set of linear filters and a point nonlinearity and are conditionally Poisson spiking. They have desirable statistical properties for fitting and have been widely used to analyze spike trains from electrophysiological recordings. However, the dynamical repertoire of GLMs has not been systematically compared to that of real neurons. Here we show that GLMs can reproduce a comprehensive suite of canonical neural response behaviors, including tonic and phasic spiking, bursting, spike rate adaptation, type I and type II excitation, and two forms of bistability. GLMs can also capture stimulus-dependent changes in spike timing precision and reliability that mimic those observed in real neurons, and can exhibit varying degrees of stochasticity, from virtually deterministic responses to greater-than-Poisson variability. These results show that Poisson GLMs can exhibit a wide range of dynamic spiking behaviors found in real neurons, making them well suited for qualitative dynamical as well as quantitative statistical studies of single-neuron and population response properties.

## 1 Introduction

Understanding the dynamical and computational properties of neurons is a fundamental challenge in cellular and systems neuroscience. A wide variety of single-neuron models have been proposed to account for neural response properties. These models can be arranged along a complexity axis ranging from detailed, interpretable, biophysically accurate models to simple, tractable, reduced functional models. Detailed Hodgkin-Huxley-style models, which sit at one end of this continuum, provide a biophysically detailed account of the conductances, currents, and channel kinetics governing neural response properties (Hodgkin & Huxley, 1952). These models can account for the vast dynamical repertoire of real neurons, but they are often unwieldy for theoretical analyses of neural coding and computation. This motivates the need for simplified models of neural spike responses that are tractable enough for mathematical, computational, and statistical analyses.

A variety of simplified dynamical models have been proposed to serve the need for mathematically tractable models, including the integrate-and-fire, Fitzhugh-Nagumo, Morris-Lecar, and Izhikevich models (FitzHugh, 1961; Nagumo, Arimoto, & Yoshizawa, 1962; Morris & Lecar, 1981; Izhikevich, 2003; Brette & Gerstner, 2005). Generally these models aim to reduce the biophysically detailed descriptions of realistic neurons to systems of differential equations with fewer variables or simplified dynamics. The one-dimensional integrate-and-fire model is arguably the simplest of these, and the simplest to analyze mathematically, but it fails to capture many of the response properties of real neurons. The two-dimensional Izhikevich model, by contrast, was specifically formulated to retain the rich dynamical repertoire of more complex, biophysically realistic models (Izhikevich, 2004).

An alternative to a mathematical notion of simplicity is the statistical property of being tractable for fitting from intracellular or extracellular physiological recordings. One well-known statistical model that satisfies this desideratum is the recurrent linear-nonlinear Poisson model, commonly referred to in the neuroscience literature as the generalized linear model (GLM) (Truccolo, Eden, Fellows, Donoghue, & Brown, 2005; Pillow et al., 2008). GLMs are closely related to generalized integrate-and-fire models such as the spike-response model, which has linear dynamics but incorporates spike-dependent feedback to capture the nonlinear effects of spiking on neural membrane potential and subsequent spike generation (Gerstner, 2001; Keat, Reinagel, Reid, & Meister, 2001; Jolivet, Lewis, & Gerstner, 2003; Pillow, Paninski, Uzzell, Simoncelli, & Chichilnisky, 2005; Gerstner, van Hemmen, & Cowan, 1996). In fact, a variant of the spike response model that incorporates noise into the spike threshold is mathematically equivalent to the models we study here (Gerstner & van Hemmen, 1992; Gerstner, 1995; Jolivet, Rauch, Lüscher, & Gerstner, 2006; Gerstner, Kistler, Naud, & Paninski, 2014). GLMs are popular for characterizing neural responses in reverse-correlation or white-noise experiments due to the tractability of likelihood-based fitting methods. Recent work has shown that GLMs can capture the detailed statistics of spiking in single and multineuron recordings from a variety of brain areas (Pillow et al., 2008; Babadi, Casti, Xiao, Kaplan, & Paninski, 2010; Calabrese, Schumacher, Schneider, Paninski, & Woolley, 2011; Weber, Machens, & Borst, 2012; Mensi et al., 2012; Pozzorini, Naud, Mensi, & Gerstner, 2013).

While several studies have shown that GLMs can successfully recapitulate various response properties of biological or simulated neurons, here we provide a more systematic study of the dynamical repertoire of the GLM. We study this issue by fitting GLMs to data from simulated neurons exhibiting a number of complex response properties. We show that GLMs can reproduce a remarkably rich set of dynamical behaviors, including tonic and phasic spiking, bursting, spike rate adaptation, type I and type II excitation, and two different forms of bistability. Furthermore, GLMs can exhibit stimulus-dependent degrees of spike timing precision and reliability (Mainen & Sejnowski, 1995) and mimic a recently reported form of greater-than-Poisson variability (Goris, Movshon, & Simoncelli, 2014).

## 2 Models of Dynamical Behaviors

### 2.1 Izhikevich Model

neuron type . | a
. | b
. | c
. | d
. | I
. | dt (ms)
. |
---|---|---|---|---|---|---|

Tonic spiking | 0.02 | 0.2 | 65 | 6 | 14 | 0.1 |

Phasic spiking | 0.02 | 0.25 | 65 | 6 | 0.5 | 0.1 |

Tonic bursting | 0.02 | 0.2 | 50 | 2 | 10 | 0.1 |

Phasic bursting | 0.02 | 0.25 | 55 | 0.05 | 0.6 | 0.1 |

Mixed mode | 0.02 | 0.2 | 55 | 4 | 10 | 0.1 |

Spike frequency adaptation | 0.01 | 0.2 | 65 | 5 | 20 | 0.1 |

Type I | 0.02 | 0.1 | 55 | 6 | 25 | 0.01 |

Type II | 0.2 | 0.26 | 65 | 0 | 0.5 | 0.01 |

Spike latency | 0.02 | 0.2 | 65 | 6 | 3.49 | 0.1 |

Resonator | 0.1 | 0.26 | 60 | 1 | 0.3 | 0.5 |

Integrator | 0.02 | 0.1 | 6 | 27.4 | 0.5 | |

Rebound spike | 0.03 | 0.25 | 60 | 4 | 5 | 0.1 |

Rebound burst | 0.03 | 0.25 | 52 | 0 | 5 | 0.1 |

Threshold variability | 0.03 | 0.25 | 60 | 4 | 2.3 | 1 |

Bistability I | 1 | 1.5 | 60 | 0 | 26.1 | 0.05 |

Bistability II | 1 | 1.5 | 60 | 0 | 26.1 | 0.05 |

neuron type . | a
. | b
. | c
. | d
. | I
. | dt (ms)
. |
---|---|---|---|---|---|---|

Tonic spiking | 0.02 | 0.2 | 65 | 6 | 14 | 0.1 |

Phasic spiking | 0.02 | 0.25 | 65 | 6 | 0.5 | 0.1 |

Tonic bursting | 0.02 | 0.2 | 50 | 2 | 10 | 0.1 |

Phasic bursting | 0.02 | 0.25 | 55 | 0.05 | 0.6 | 0.1 |

Mixed mode | 0.02 | 0.2 | 55 | 4 | 10 | 0.1 |

Spike frequency adaptation | 0.01 | 0.2 | 65 | 5 | 20 | 0.1 |

Type I | 0.02 | 0.1 | 55 | 6 | 25 | 0.01 |

Type II | 0.2 | 0.26 | 65 | 0 | 0.5 | 0.01 |

Spike latency | 0.02 | 0.2 | 65 | 6 | 3.49 | 0.1 |

Resonator | 0.1 | 0.26 | 60 | 1 | 0.3 | 0.5 |

Integrator | 0.02 | 0.1 | 6 | 27.4 | 0.5 | |

Rebound spike | 0.03 | 0.25 | 60 | 4 | 5 | 0.1 |

Rebound burst | 0.03 | 0.25 | 52 | 0 | 5 | 0.1 |

Threshold variability | 0.03 | 0.25 | 60 | 4 | 2.3 | 1 |

Bistability I | 1 | 1.5 | 60 | 0 | 26.1 | 0.05 |

Bistability II | 1 | 1.5 | 60 | 0 | 26.1 | 0.05 |

### 2.2 Generalized Linear Model

The GLM is a regression model typically used to characterize the relationship between external or internal covariates and a set of recorded spike trains. In systems neuroscience, the label “GLM” often refers to an autoregressive point process model, a model in which linear functions of stimulus and spike history are nonlinearly transformed to produce the spike rate or conditional intensity of a Poisson process (Truccolo et al., 2005; Pillow et al., 2008).

Unlike classical deterministic models like Hodgkin-Huxley and integrate-and-fire, the GLM is fundamentally stochastic due to the assumption of conditionally Poisson spiking. However, this stochasticity is helpful for fitting purposes because it assigns graded probabilities to firing events and allows for likelihood-based methods for parameter fitting (Paninski, Pillow, & Simoncelli, 2004; Pillow et al., 2005). In fact, the Poisson GLM comes with a well-known guarantee that the log-likelihood function is concave for suitable choices of nonlinearity (Paninski, 2004). This means we can be assured of approaching a global optimum of the likelihood function via gradient ascent for any set of stimuli and spike trains (barring any numerical issues that may complicate achieving the actual maximum for certain data sets; cf. Zhao & Iyengar, 2010). This guarantee does not hold for stochastic formulations of most nonlinear biophysical models, including the Izhikevich model. Moreover, despite its stochasticity, the GLM can produce highly precise and repeatable spike trains in certain parameter settings, as we will demonstrate.

## 3 GLMs Capture a Wide Array of Complex Dynamical Behaviors

We fit GLMs to data simulated from Izhikevich neurons set up to exhibit a range of different qualitative response behaviors. In the following, we describe these behaviors in detail, beginning with simpler behaviors, such as tonic spiking and bursting (which have already been demonstrated in previous work—e.g., Gerstner & van Hemmen, 1992; Jolivet et al., 2006) in order to build intuition for the GLM's basic capabilities, and then move on to more complex behaviors (such as bistability) and questions of spike timing reliability and precision.

### 3.1 Tonic Spiking

We first examined an Izhikevich neuron tuned to exhibit tonic spiking (see Figures 2A and 2B; see Table 1 for parameters). When presented with a step input current, the Izhikevich neuron responds with a few high-frequency spikes and then settles into a regular firing pattern that persists for the duration of the step (see Figure 2B). This response pattern resembles that of a deterministic Hodgkin-Huxley or integrate-and-fire neuron, albeit with an added transient burst of spikes at stimulus onset.

We simulated the Izhikevich neuron's response to a series of step currents and used the resulting training data to perform maximum-likelihood fitting of the GLM parameters . The estimated stimulus filter is biphasic (see Figure 2C), resulting in a large transient response to stimulus onset, and has a positive integral, ensuring a sustained positive response to a current step (see Figure 2E, blue trace). The estimated post-spike filter (see Figure 2D), by contrast, has a large negative lobe that provides recurrent inhibition after every spike, enforcing a strong relative refractory period. The stimulus filter and post-spike filter output (shown together for a single trial in Figure 2E) are summed together and exponentiated to obtain the conditional intensity (see Figure 2F), also known as the instantaneous spike rate.

For this stimulus, the intensity rises very quickly once decays, which occurs approximately 50 ms after the previous spike. Note that the output of the stimulus filter is identical on each trial, whereas the output of the post-spike filter varies from trial to trial because of variability in the exact timing of spikes. However, because the rising phase of the conditional intensity is so rapid, spiking is virtually certain within a small time window sitting at a fixed latency after the previous spike time. The combination of a strong excitatory drive from the stimulus filter and a strong suppressive drive from the post-spike filter produces precisely timed spikes across trials, allowing the GLM to closely match the deterministic firing pattern of the Izhikevich neuron (see Figure 2G).

### 3.2 Bursting

We next examined multispike bursting, a more complex temporal response pattern that requires dependencies beyond the most recent interspike interval. Once again, we simulated responses from an Izhikevich neuron tuned to exhibit tonic bursting (see Figures 3A and 3B) and used the resulting data to fit GLM parameters (see Figures 3C and 3D). The estimated stimulus filter is biphasic with a larger positive than negative lobe, which drives rapid spiking at stimulus onset and generates sustained drive during an elevated stimulus (see Figure 3E, blue trace). The post-spike filter has an immediate negative component that creates a relative refractory period after each spike and an even more negative mode after a latency of 40 ms; the accumulation of these negative components over multiple spikes gives rise to a sustained suppression of activity between bursts.

The GLM captures the bursting behavior of the Izhikevich neuron with high precision, including the fact that the first burst after stimulus onset contains a different pattern of spikes than subsequent bursts do. This difference arises from precise interactions between the stimulus and post-spike filter outputs. During the first burst, fast spiking arises from an interplay between monotonically increasing stimulus filter output (see Figure 3E, blue) and tonic decrements induced by the post-spike filter after each spike (see Figure 3E, red). After each spike, the post-spike filter reduces the conditional intensity by a fixed decrement, but this decrement is soon overwhelmed by the rising wave of input from the stimulus filter, which creates a rapid rise and fall of the conditional intensity time-locked to each Izhikevich neuron spike time. The pattern continues until accumulated contributions from the delayed negative lobes of post-spike filter overwhelm those from the stimulus filter and the burst terminates. Subsequent bursts are governed by a somewhat different interplay between stimulus and post-spike filter outputs: bursts sit on a rising phase of the conditional intensity due to the removal of suppression from the previous burst. This rise is more gradual than the drive induced by stimulus onset and results in bursts with longer interspike intervals and fewer spikes per burst, but the resulting spike pattern is nonetheless captured with high precision and reliability from trial to trial.

### 3.3 Bistability

Bistability refers to the phenomenon in which there are multiple stable response modes for a single input condition. A common form of bistability observed in real neurons is the ability to inhabit either a tonically active state or a silent state for a given level of current injection. The Izhikevich model can exhibit this form of bistability, wherein a brief positive current pulse is sufficient to kick it between states: the neuron can inhabit a silent state in the absence of stimulation, but a brief positive current pulse kicks it into a tonically active state, and an appropriately timed positive pulse kicks it back to the silent state (see Figures 4A and 4B).

We fit a GLM to spike trains simulated from such a bistable Izhikevich neuron and found that the fitted GLM can capture bistable behavior of the original model with high accuracy (see Figures 4C to 4G). When stimulated during the silent state, the GLM emits a spike due to the positive output of the stimulus filter (see Figure 4E, blue), and tonic firing ensues due to a positive lobe in the post-spike filter that causes self-excitation at a fixed latency of approximately 10 ms after the previous spike (see Figure 4E). The GLM returns from the active state to the silent state when a positive stimulus pulse synchronizes negative lobes of the stimulus and post-spike filters. Only the combination of these two negative drives is strong enough to shut off spiking. Without the negative drive created by previous spikes (appearing in the post-spike filter at a latency of 15 ms after a spike), suppression from the negative lobe of the stimulus filter is not strong enough to prevent spiking. As with the tonic spiking neuron, the interaction between the stimulus and post-spike filters generates rapid rises in the conditional intensity (see Figure 4F), leading to precisely timed spikes that mimic those of the deterministic Izhikevich neuron (see Figure 4G).

This Izhikevich neuron (with the same parameters) can also exhibit a second form of bistability in which the return to the silent state from the active state is induced by a negative instead of a positive current pulse. We performed a similar fitting exercise and found that the GLM is also able to reproduce this behavior. Firing is initiated and maintained by a similar mechanism as the first form of bistability, but firing offset occurs due to the fact that a negative stimulus pulse creates immediate negative output from the stimulus filter, which suppresses firing during the time when a spike would have occurred due to spike history filter input. Tonic firing is extinguished more rapidly in this second form of bistability than the first (see Figure 6).

### 3.4 Type I and Type II Firing

Neurons have been classified as exhibiting either type I or type II dynamics based on the shape of their firing rate versus intensity (F-I) curve. Type I neurons can fire at arbitrarily low rates for low levels of injected current, whereas type II neurons have discontinuous F-I curves that arise from an abrupt transition from silence to a finite nonzero firing rate as the level of injected current increases (Hodgkin, 1948). We simulated Izhikevich neurons that exhibit each of these response types using published parameters. Inputs consisted of 500 ms current steps of varying amplitude. The resulting F-I curves for the Izhikevich type I and type II neurons are shown in black in Figures 5A and 5B, respectively. We fit GLMs using data from each Izhikevich neuron and found that the fitted GLMs capture the two response types with high temporal precision. The corresponding F-I curves are shown in gray in Figure 5 and accurately mimic the behaviors of the Izhikevich neuron. Similar F-I curves have been demonstrated previously in Gerstner and van Hemmen (1992) and Mensi et al. (2012).

The only discrepancy between the Izhikevich and GLM neurons occurs for the type II cell at input amplitudes near the Izhikevich neuron's threshold. On some trials when the input amplitude falls below this threshold, the GLM jumps into a tonic firing state for the duration of the stimulus. Similarly, on some trials when the input amplitude falls above this threshold, the GLM fails to initiate firing. This is unsurprising given the stochastic nature of the GLM. Importantly, the GLM never fires at a low rate, but rather abruptly transitions from no firing to firing at a baseline level of 25 Hz, reflecting type II behavior.

### 3.5 Additional Behaviors

We fit GLMs to every dynamical behavior considered in Izhikevich (2004) with the exception of purely subthreshold behaviors, since GLM fitting uses spike trains and does not consider subthreshold responses. The full suite of behaviors is shown in Figure 6, with responses of the Izhikevich neurons in black and spike responses of the GLM in gray. This list includes tonic and phasic spiking, tonic and phasic bursting, mixed mode firing, spike frequency adaptation, type I and type II excitability, two different forms of bistability, and several others that depend primarily on the shape of stimulus filter. Several additional behaviors that can be captured by a GLM are not depicted in Figure 6 as they can be achieved by a trivial manipulation of the stimulus filter; for example, inhibition-induced bursting can be achieved by simply flipping the sign of the stimulus filter for the bursting neuron shown in Figure 6C. Previous work has shown the Izhikevich neuron to be capable of producing 18 distinct spiking behaviors (Izhikevich, 2004), and we found that all can also be produced by a GLM.

### 3.6 Systematic Variation of Filter Amplitudes

We next considered what happens to the behaviors produced by a GLM as some aspect of the filters is systematically varied. To do so, we created stimulus and post-spike filters composed by linear combinations of two basis filters and then systematically varied the amplitude of one basis filter while holding the other fixed. (See the appendix for details.) Figure 7 shows the phase space of qualitative spiking behaviors obtained at different points in this 2D filter space.

When the stimulus filter has a strongly negative component (center panel, bottom), a positive stimulus pulse does not produce enough driving force to cause the neuron to spike at all (quiescent). As the amplitude of this component of the stimulus filter is increased, the neuron receives stronger and stronger input and is driven first to spike once or twice (phasic spiking) and eventually to emit a burst of spikes (phasic bursting). The stimulus filter largely drives changes between these behaviors, with the additional detail that a strongly negative post-spike filter component is able to inhibit a burst that would otherwise occur (upper left corner of “phasic spiking” region). As the stimulus filter component becomes still more positive and the stimulus filter transitions from being biphasic to more monophasic, it produces a positive driving force for the duration of the stimulus step rather than just the onset. This causes the neuron to fire for the duration of the step.

Importantly, the post-spike filter here determines the nature of this sustained firing. A post-spike filter that is purely negative beginning at short timescales (top, left side) mimics a relative refractory period, inhibiting additional spikes for a short window following each elicited spike and resulting in tonic spiking. If the post-spike filter is only weakly negative at short timescales while being more strongly negative at longer timescales, this creates multiple timescales in the neuron's response (top, right side). At short timescales, there is little inhibition from each spike (beyond the absolute refractory period), so additional spikes may occur. Over longer timescales, inhibition is accumulated over multiple spikes, which eventually shuts off spiking. After a brief window of no spikes, the inhibition is relaxed and spiking begins again until enough inhibition is accumulated to shut spiking off. This cycle results in tonic bursting for the duration of the step. In the extreme case where the post-spike filter is actually biphasic (top, far right), each spike promotes additional spikes on short timescales, leading to highly regular timing of spikes within bursts (see Figure 7C). This set of behaviors could be achieved by simply sweeping over the amplitude of a single basis vector in each filter. Incorporating shifts to the basis vectors or additional basis vectors would likely be necessary to achieve more complex behaviors, such as bistability.

Although we have drawn clear borders at the transition between behaviors, these transitions in fact occur gradually. Near the border between phasic spiking and phasic bursting, for example, there will be some trials where a single spike is produced and other trials where a burst is elicited. We have indicated the behavior that is produced most frequently here for simplicity. The transition from tonic bursting to tonic spiking also occurs gradually, with the near perfectly regular bursting breaking down into more irregular firing until no apparent bursts are produced. If the post-spike filter is made even more negative than the range explored in this figure, the timing of tonic spiking becomes nearly perfectly regular as well. This is easily explained by the fact that as the post-spike filter becomes more and more negative, it imposes stronger refractoriness on the cell, which results in more regular spike timing. As the post-spike filter component amplitude is changed, there is a gradual change from precisely timed bursts, to irregular firing, to precisely timed tonic spiking. In sections 4 and 5, we further explore questions of spike timing precision in the GLM.

### 3.7 Generalization to New Stimuli

The behaviors shown in Figure 6 were all fit using stimuli that probe only a small range of the possible behaviors of the neuron. For example, many were probed using only a single step height. A natural question that arises is therefore: how well will these fitted GLMs generalize to predict the Izhikevich neuron responses to new stimuli? We examined this question for three canonical Izhikevich neurons from our study (see Figure 8). We first generated responses from a regular spiking Izhikevich neuron using three step heights and one noise stimulus (see Figures 8A and 8B, top). We then simulated responses to these stimuli using the GLM fit only to the intermediate step height (see Figure 8B, middle). (This is the same fit as Figure 2.) The GLM responses nearly perfectly capture the Izhikevich responses for the original stimulus, and the GLM maintains regular firing patterns for the other step heights. However, the GLM's firing rate is too high for the small step and too low for the large step. While the GLM accurately captures some firing events for the noise stimulus, the firing rate is overall too high. We next refit a GLM to responses from all three step heights using the same set of basis vectors as the original fit (see Figure 8, bottom). The responses to all three step heights are captured nearly perfectly. Additionally, the response to the noise stimulus is much more similar to that of the Izhikevich neuron. Although there is one firing event that occurs at a delay, the regular spiking GLM fit on an enriched stimulus set is better able to generalize to the noise stimulus.

We performed this same test for neurons showing phasic bursting (see Figure 8C) and tonic bursting (see Figure 8D). (Note that although some responses of the phasic bursting Izhikevich neuron are not actually phasic, we retain the naming convention given to this set of parameters in the original paper.) For phasic bursting, the GLM fit on additional step sizes improves the accuracy of responses to the smallest and largest steps while decreasing the accuracy of responses for the original step of intermediate size. There is marked improvement in the accuracy of responses to the noise stimulus, with many firing events being accurately captured and the GLM no longer exhibiting runaway excitation. For tonic bursting, the refit GLM retains bursting behavior but fails to even capture responses for the steps on which it was trained. As noted above, when refitting, we used the same set of basis vectors as the initial fits for fair comparison. It is possible that by increasing the number of basis vectors used or tuning their properties that better fits to all stimuli might be achieved.

Taken together, these results show that while GLMs might retain some characteristic response features (such as bursting) when probed with new stimuli, they often have limited ability to generalize beyond stimuli on which they are directly fit.

## 4 Spiking Precision and Reliability

A noteworthy feature of the spike trains of the GLM neurons considered above is their high degree of spike timing precision and reliability across trials. This precision arises from the fact that the conditional intensity (or instantaneous spike rate) rises abruptly at spike times (due to filter outputs passing through a rapidly accelerating exponential nonlinearity) and decreases immediately after each spike due to suppressive effects of the post-spike filter. By contrast, a Poisson GLM without recurrent feedback, more commonly known as a linear-nonlinear Poisson (LNP) cascade model, cannot produce temporally precise spike responses to a constant stimulus because its output is constrained to be a Poisson process.

Real neurons, however, seem to be capable of both response modes: they emit precisely timed spikes in some settings and highly variable spike trains in others. A seminal paper by Mainen and Sejnowski (1995) illustrated this duality by showing that spike responses to a constant current exhibit substantial trial-to-trial variability, whereas responses to a rapidly fluctuating injected current are precise and repeatable across trials. Deterministic models like the Izhikevich model cannot, of course, mimic this property because their spikes are perfectly reproducible for any stimulus. (A stochastic version of the Izhikevich model with an appropriate level of injected noise could likely overcome this shortcoming, however. See Schneidman, Freedman, and Segev (1998) for a similar case in a Hodgkin-Huxley neuron.) Here we show that the GLM naturally reproduces the same form of stimulus-dependent changes in precision and reliability observed in real neurons. Figure 9 shows that a single GLM (with parameters identical to those fit to the tonic-spiking Izhikevich neuron, shown in Figure 2) produces irregular spiking in response to a constant stimulus with low-to-intermediate amplitude and precisely timed, reliable spikes in response to a stimulus with large, rapid fluctuations.

## 5 GLMs Can Produce Super-Poisson Variability

We have shown that GLM neurons can reproduce the high degree of spike timing precision found in real neurons stimulated with injected currents. However, a variety of studies have reported that neurons exhibit overdispersed responses, or greater-than-Poisson spike count variability in response to repeated presentations of a sensory stimulus (Tomko & Crapper, 1974; Tolhurst, Movshon, & Dean, 1983; Shadlen & Newsome, 1998; Goris et al., 2014). A prominent recent study from Goris et al. (2014) showed that the degree of overdispersion grows with mean spike count, so that the Fano factor (variance-to-mean ratio) is an increasing function of spike rate. They proposed a doubly stochastic model to account for this phenomenon, in which the rate of a Poisson process is modulated by a slowly fluctuating stochastic gain variable . For each trial, is drawn from a gamma distribution with mean 1 and variance . (See the appendix.)

We sought to determine if a GLM with spike-history dependence can also account for the mean-dependent overdispersion found in neural responses. To test this possibility, we simulated spike trains from the doubly stochastic model of Goris et al. (2014) for three different settings of the overdispersion factor with the same mean spike rate (100 spikes/s), and fit a GLM to the spike trains associated with each value of . We then simulated responses from each GLM to 500 ms pulses at a number of different input intensities, with each point in Figure 10 corresponding to a different intensity.

We found that the GLM can indeed match the qualitative behavior of the Goris et al. (2014) model, giving approximately Poisson responses at low spike rates and increasingly overdispersed responses at higher rates (see Figure 10C). To match the data from larger values of overdispersion factor (see the darker curves in Figure 10C), the GLM relies on increasing amounts of self-excitation from the post-spike filter but exhibits no changes in stimulus filter (see Figures 10A and 10B). The filters here do not include an absolute refractory period, as the original model does not incorporate one. However, similar results can be achieved when a refractory period is enforced in the training data. This will result in post-spike filters with strongly negative lobes on short timescales, which impose refractoriness but otherwise similar filters to those in Figure 10. Unlike models with purely suppressive spike history effects, which capture effects due to refractoriness and reduce variability in firing (e.g., Berry and Meister, 1998; Uzzell, 2004), here we show that allowing spike history effects to be excitatory can result in increased variability. While the former might be suitable for early sensory areas such as the retina, the latter better captures the super-Poisson variability observed in higher visual areas.

Intuitively, the GLM generates overdispersed spike counts because of dependencies the spike-history filter induces between early and late spikes during a trial. If the GLM neuron generates a larger-than-average number of spikes early in a trial, the positive post-spike filter produces a higher conditional intensity (and hence more spiking) later in the trial; conversely, if a neuron emits fewer-than-average spikes early in a trial, the conditional intensity will be lower later in the trial (yielding less spiking). The Goris et al. (2014) model can be seen to capture similar dependencies between early and late spikes via the stochastic gain variable , which is constant during a trial but independent across trials. Thus, it is reasonable to view in the Goris et al. (2014) model as a proxy for the accumulated self-excitation from spike-history filter outputs under a GLM. We note, however, that attempts to drive the GLM to higher levels of overdispersion (e.g., Fano factors significantly larger than 3) often resulted in runaway self-excitation, indicating that GLMs may require additional mechanisms to maintain stability in order to produce highly overdispersed responses through recurrent excitation alone (Gerhard, Deger, & Truccolo, 2017; Hocker & Park, 2017). An alternative mechanism for generating overdispersed responses with GLMs is through the addition of latent stochastic inputs (Rabinowitz, Goris, Cohen, & Simoncelli, 2015), an avenue we have not explored here.

## 6 Discussion

We have shown that recurrent Poisson GLMs can capture an extensive set of behaviors exhibited by biological neurons, including tonic and phasic spiking, bursting, spike frequency adaptation, type I and type II behavior, and bistability. GLMs can also reproduce widely varying levels of response stochasticity, ranging from precisely timed spikes with negligible trial-to-trial variability, to substantially super-Poisson spike count variability. We have also shown that like real neurons, GLMs can exhibit irregular firing in response to a constant stimulus but precise and repeatable firing patterns in response to a temporally varying stimulus. Thus, generalized linear models are able to capture a rich array of spiking behaviors like many dynamical models while remaining tractable to fit to neural data.

### 6.1 Relationship to Previous Work

GLMs have strong connections to a number of other models. It is particularly worth noting the connection between GLMs and generalized integrate-and-fire models, such as the spike response model (SRM) extensively studied by Gerstner and colleagues (Gerstner & van Hemmen, 1992; Gerstner, 2001). These models draw on much earlier work that incorporated a variable threshold that depends on spiking history (Weiss, 1966; Geisler & Goldberg, 1966). The SRM includes a membrane filter (analogous to the stimulus filter here) and both a spike afterpotential and moving threshold (which can be combined and are analogous to the spike history filter here). In its simplest formulation, the SRM is a deterministic model. Although the threshold for spiking can shift as a function of spike history, a spike will occur precisely at each threshold crossing.

Extensions of the model have incorporated so-called ``escape noise,'' where spiking no longer occurs deterministically at threshold crossings, but rather the probability of spiking depends on the distance of the membrane voltage to threshold (Gerstner & van Hemmen, 1992; Jolivet et al., 2006). This variant of the SRM is in fact a GLM, and it is therefore worthwhile to consider how previous work investigating the SRM with escape noise relates to our results here. Early work demonstrated that such a model was capable of producing responses with high temporal precision, including tonic spiking and tonic bursting (Gerstner & van Hemmen, 1992), though demonstrating the range of behaviors that could be produced by the SRM was not the focus of this study. Additional work demonstrated that the model could produce highly repeatable spike trains to a noisy stimulus (similar to Figure 9B, though no comparison of irregular firing in response to a constant stimulus was shown; Jolivet et al., 2006). Other studies have shown that the SRM can capture the detailed statistics of neural responses (Mensi et al., 2012; Pozzorini et al., 2013). Further, many of these studies show that the spike response model can be used to capture not only the relationship between an external stimulus and a neuron's response but also to faithfully capture the relationship between intracellularly recorded neural responses and injected current (Jolivet et al., 2006; Pozzorini et al., 2013).

### 6.2 Limitations

Despite their many advantages, GLMs have several limitations that bear further discussion. First, they often do not generalize well across stimulus distributions; models fit with a particular set of stimuli often do not accurately predict responses to stimuli with markedly different statistics (e.g., stimuli with large changes in mean or variance, or white noise versus naturalistic stimuli) (Heitman et al., 2016).

Second, GLMs often lack clear interpretability in terms of underlying mechanisms. This stands in contrast to dynamical models designed to capture specific biophysical variables and processes. In the two-dimensional Izhikevich model, for example, one variable () represents the neuron's membrane potential, and the other variable () can be understood as a membrane recovery variable, which reflects K^{+} channel activation and Na^{+} channel inactivation. Despite the fact that the GLM filters do not represent specific biophysical variables, in some cases they can still provide insight into underlying biological processes. For example, recent work provides an interpretation of the GLM as a synaptic conductance-based model with linear subthreshold dynamics (Latimer, Chichilnisky, Rieke, & Pillow, 2014). Work on the SRM has shown that by dividing the effects of spike history into a dynamic threshold and a spike afterpotential, one can in fact measure their separate contributions with intracellular recordings of a neuron's subthreshold voltage; the spike afterpotential can be observed directly in this voltage trace, while the effects on threshold can be estimated indirectly by noting the absence of firing (Mensi et al., 2012; Pozzorini et al., 2013; Gerstner et al., 2014).

A third known limitation of GLMs is that they lack the flexibility to capture some nonlinear response properties of real spike trains. For example, as point neuron models, GLMs do not reflect the fact that neurons often receive spatially segregated inputs on the dendritic tree, and these inputs can be processed separately and combined nonlinearly (London & Häusser, 2005). Some extensions of the GLM that incorporate nonlinear inputs and multiple subunits (Fitzgerald, Rowekamp, Sincich, & Sharpee, 2011; Rajan, Marre, & Tkačik, 2013; Park, Archer, Priebe, & Pillow, 2013; Cui, Liu, Khawaja, Pack, & Butts, 2013; McFarland, Cui, & Butts, 2013; Ahrens, Linden, & Sahani, 2008; Williamson, Sahani, & Pillow, 2015) may begin to address this issue but certainly fall short of capturing the full complexity of dendritic processing. For the range of dynamical behaviors considered here, however, we did not find these extensions to be necessary.

For all results shown, we used a GLM with an exponential nonlinearity. To test the dependence of our results on the form of the nonlinearity, we also fit GLMs to several of the behaviors with a “soft-rectifying” nonlinearity given by . This function grows only linearly for large input values, but still has an exponential decay on its left tail and remains in the family of nonlinearities (convex and log-concave) for which the GLM log-likelihood is provably concave (Paninski, 2004). For the behaviors tested (tonic spiking, tonic bursting, phasic spiking, and phasic bursting), our results were similar to those with an exponential nonlinearity, though generally not as temporally precise. This increased precision is likely due to the fact that an exponential nonlinearity rises more steeply than a linear-rectifying function, causing the conditional intensity to accelerate more rapidly from a low probability to a high probability of spiking regime. Past studies have found that responses of both retinal ganglion cells and neocortical pyramidal neurons are well described by a GLM with exponential nonlinearity (Pillow et al., 2008; Jolivet et al., 2006).

It is worth noting that for many of the dynamic behaviors studied here, the GLM parameters were not strongly constrained by the training data. (See section A.4 in the appendix.) Slight changes in the model parameters did not produce noticeable changes in response, at least for the stereotyped range of input currents and output spike patterns considered. The filter parameters were therefore only weakly identifiable, which corresponds to a likelihood function with a very gradual falloff along certain directions in parameter space. This uncertainty potentially complicates interpretation of the filters in terms of functions performed by the underlying biophysical mechanism. Conversely, it reveals that the suite of behaviors that Izhikevich and others considered can be achieved by a range of different GLMs and that a richer set of input-output patterns is needed to identify a unique set of GLM parameters.

## 7 Conclusion

The GLM has the ability to mimic a wide range of biophysically realistic behaviors exhibited by real neurons. Although it is clear there are some forms of nonlinear behavior it cannot produce, such as frequency-doubled responses of cat Y cells or V1 complex cells to a contrast-reversing grating, our work provides an existence proof for its ability to exhibit an important range of response types considered previously only in biophysics and applied math modeling literature. Moreover, by considering response stochasticity as another dimension along which real neurons vary, we have shown that that the GLM can generate response characteristics ranging from quasi-deterministic to greater-than-Poisson variability. The GLM therefore provides a flexible yet powerful tool for studying the dynamics of real neurons and the computations they carry out.

## Appendix

Matlab code used to generate example responses from Izhikevich neurons and to fit GLMs to these responses is available in a Github repository (https://github.com/aiweber/GLM_and_Izhikevich).

### A.1 Izhikevich Model Simulations

To generate training data for fitting the GLM, we simulated responses from an Izhikevich model (Izhikevich, 2003) with parameters set to published values given for each behavior in Izhikevich (2004) (see Table 1; parameter values can be found at http://www.izhikevich.org/publications/izhikevich.m). For each behavior, we generated approximately 20 seconds of training data using the forward Euler method with fixed time step size () given in Table 1. It should be noted that in some cases, published parameter values did not produce the desired qualitative behavior. In these cases, we tuned the simulation parameters to achieve the desired behavior. Parameters marked with an asterisk in Table 1 indicate those that differ from published values for the corresponding behavior in Izhikevich (2004). Additionally, some behaviors of the Izhikevich neuron are not robust to small changes in stimulus timing, stimulus amplitude, or time step of integration. In particular, we found bistability to be highly dependent on the precise stimulus timing (onset and duration), stimulus amplitude, and integration window. We tuned these values by hand to produce the desired behavior.

### A.2 GLM Fitting and Simulations

### A.3 Systematic Variation of Filter Amplitudes

In order to more carefully examine the transitions between different behaviors as the stimulus and post-spike filter change, we systematically varied the amplitude of individual filter components and observed the behavior produced. Each filter was parameterized with two components. The amplitude of one was fixed, while the amplitude of the other was varied. For the stimulus filter, the amplitude of the second component was varied (1.5 to 0.25). This allowed us to transition from monophasic to biphasic filters. The amplitude of the first component was set to be positive (1), creating an ON filter appropriate for a positive step stimulus. For the post-spike filter, the amplitude of the first component was varied (1 to 1). The amplitude of the second component was set to be negative (3), ensuring that spiking would be suppressed on longer timescales. For the post-spike filter, we also imposed an absolute refractory period of 5 ms. Finally, we included a negative baseline drive ( = 1) to suppress spontaneous spiking so that the baseline firing rate was zero.

We simulated responses to 25 identical step stimuli for each set of filters and then classified the behaviors as quiescent, phasic spiking, phasic bursting, tonic spiking, or tonic bursting. The most commonly observed behavior over the 25 repetitions is depicted in Figure 7. Responses were classified in the following way. If no spikes were elicited in the first 200 ms of stimulus presentation and fewer than five spikes were elicited during the final 10 seconds of stimulus presentation, the behavior was classified as quiescent. If at least one spike was elicited in the first 200 ms following stimulus onset and fewer than five spikes were elicited during the final 10 seconds of stimulus presentation, the behavior was classified as phasic. Phasic firing patterns were further classified into phasic spiking if only one or two spikes were elicited in the first 200 ms, and phasic bursting if three or more spikes were elicited in the first 200 ms. The remaining responses were classified as either tonic spiking or tonic bursting in the following manner. Interspike interval distributions were fit with a gaussian mixture distribution using Matlab's gmdistribution function. We fit both a single gaussian distribution and a mixture of two gaussians and then compared the Akaike information criterion (AIC) values to determine whether the ISI distribution was better fit as a unimodal distribution or a bimodal distribution, with a lower AIC indicating better fit. If , the spike train was classified as tonic spiking; otherwise, it was classified as tonic bursting. (We added the 0.9 factor to create a more stringent standard for what is classified as bursting activity so that the responses that fall into this category are strongly bimodal distributions that would be readily identified as bursting. Slightly altering the value of this factor or eliminating it entirely gives the same qualitative results, but merely shifts the boundary in Figure 7 between the tonic spiking and tonic bursting regions.)

### A.4 Sensitivity to Changes in Parameter Values

We wished to investigate how well constrained different features of our fit GLMs were. To do so, we calculated the eigendecomposition of the Hessian matrix of the likelihood function: . The Hessian matrix provides a local quadratic approximation to the likelihood function, with eigenvectors pointing along the principal axes and length of these axes proportional to . Thus, larger-magnitude eigenvalues indicate greater curvature (i.e., shorter axes) and better-constrained directions, while smaller-magnitude eigenvalues indicate lower curvature and more poorly constrained directions. Results of this analysis are shown in Figure 11 for both a regular spiking neuron (left) and tonic bursting neuron (right). For both neurons, the least constrained directions correspond to eigenvectors similar in shape to the best-fit filters or the absolute refractory period. As such, perturbations in these directions do not result in large changes in the behavior. Perturbations along eigenvectors corresponding to the most constrained directions would result in significant changes to the filter shape. The difference in scale between panel A and panel D indicates that overall, parameters for the tonic bursting neuron are more constrained than those for the regular spiking neuron.

### A.5 Doubly Stochastic Model with Super-Poisson Variability

For the results shown in Figure 10, we simulated data from the negative binomial distribution with a single mean rate (100 spikes/s) at three different gain variances . Spike counts were drawn i.i.d. across trials, with spike times distributed uniformly within each trial to generate spike trains suitable for GLM fitting. We used these spike trains to fit a GLM to the data associated with each value of , with an assumed constant input current for each trial.

## Acknowledgments

We thank Adrienne Fairhall, Sara Solla, and James Fitzgerald for helpful comments and discussions.