Biol. Cybern. 77, 289–295 (1997)

c

� Springer-Verlag 1997

**Satish Iyengar ^{}**1

Received: 16 September 1996 / Accepted in revised form: 2 July 1997

**Abstract. **Spike trains from neurons are often used to make inferences about the underlying processes that generate the spikes. Random walks or diffusions are commonly used to model these processes; in such models, a spike corresponds to the ﬁrst passage of the diffusion to a boundary, or ﬁring threshold. An important ﬁrst step in such a study is to ﬁt families of densities to the trains’ interspike interval histograms; the estimated parameters, and the families’ goodness of ﬁt can then provide information about the process leading to the spikes. In this paper, we propose the generalized inverse Gaussian family because its members arise as ﬁrst passage time distributions of certain diffusions to a constant boundary. We provide some theoretical support for the use of these diffusions in neural ﬁring models. We compare this family with the lognormal family, using spike trains from retinal ganglion cells of goldﬁsh, and simulations from an integrate-and-ﬁre and a dynamical model for generating spikes. We show that the generalized inverse Gaussian family is closer to the true model in all these cases.

The analysis of neural spike trains has a long history (Tuckwell 1988). The study of spike trains has several purposes: for example, to draw inferences about the nature of the inputs to the neuron, and the processes within the neuron that convert those inputs into impulses. There are two main approaches for such analyses. The ﬁrst is to formulate a stochastic model for the neuron’s activity and inputs, derive the distribution of the interspike interval (ISI), and use it to ﬁt the model to data (Gerstein and Mandelbrot 1964; Inoue et al. 1995). The stochastic process is often modeled as a random walk or diffusion, and the spike corresponds to its ﬁrst passage time to a boundary, or ﬁring threshold. When such a stochastic model is not known, a second approach is to ﬁt standard families of distributions, such as the lognormal or gamma, using maximum likelihood or some ad hoc technique. In this case, inferences about the underlying

*Correspondenc*e *to*: S. Iyengar (Tel: +1 412 624 8341, Fax: +1 412 624 8397, e-mail: si@stat.pitt.edu)

process are more difﬁcult, and at times not attempted at all (Tuckwell 1988; Levine 1991).

In this paper, we study the adequacy of the generalized inverse Gaussian family for ISI data. We propose this family for several reasons. First, for certain values of the parameters of this family, its members are the ﬁrst passage time distributions of certain diffusion processes to a constant boundary. Also, it is an exponential family which contains commonly used subfamilies such as the gamma, inverse Gaussian, and their reciprocals. Using spike trains from retinal ganglion cells of goldﬁsh, we compare this family with the lognormal family, which is commonly used to ﬁt ISI histograms. For problems of this sort, it is typically not known that the true model (the one that represents the process which generates the data) belongs to one of the families under consideration. We therefore use a quasi-maximum likelihood approach to identify the family that is closer to the true model, using the Kullback-Leibler information criterion (Nishii 1988). We also apply this technique to spike trains simulated from an integrate-and-ﬁre model, and from a dynamical model.

In Sect. 2, we present the spike train data, give the details of the two families, and describe the neurophysiological signiﬁcance of the diffusion processes underlying the generalized inverse Gaussian family. We then describe the method of model comparison, and derive the test statistics needed. In Sect. 3, we give the results of the tests, and show that the generalized inverse Gaussian family is the closer model for all cases. We then conclude with a discussion of our ﬁndings in Sect. 4.

The simulated data are from two models of spike generation. The ﬁrst is Stein’s integrate-and-ﬁre model (Tuckwell 1988). We use noise uniformly distributed on the interval [0, 1] to perturb a mean bias signal to create the integrated variable. The summary statistics of the noise, i.e. [mean, standard deviation, skewness, kurtosis], for the two runs are [0.501, 0.288, 0.001, *−*1*.*21] and [0.501, 0.289, 0.001, *−*1*.*18], respectively; these values are close to their expected values. The integration proceeds at 1-ms time steps until a selected ﬁring threshold is reached. This model uses a leaky integrator with time constant τ = 5 ms, and approximately 30 000 noise inputs **Table 1. Summary statistics for integrate-and-ﬁre simulations **

Run: | 1 | 2 |
---|---|---|

Spikes: | 1159 | 1102 |

Mean: | 25.9 | 26.4 |

Median: | 22.0 | 23.0 |

SD: | 15.9 | 11.9 |

CV: | 0.61 | 0.45 |

**Table 2. Summary statistics for dynamical model simulations **

**Run: 12 **

**Spikes: 6383 6364 **

**Mean: 3.13 3.14 **

**Median: 2.89 2.91 **

**SD: 1.06 1.02 **

**CV: 0.34 0.32 **

**to produce more than 1100 simulated spikes. Table 1 contains summary statistics for these two simulated spike trains. The ISI means are around 26 ms, the ISI histograms are skewed to the right, and their coefﬁcients of variation (CVs) are 0.45 and 0.61. Figure 1 shows a smoothed histogram of the latter, using a 1-ms binwidth. **

**The second set of simulated data is from two runs of a dynamical model which is a variant of a model due to Morris and Lecar 1981; for background and a description of the simulation methods. see Ermentrout 1996; Ermentrout and Rinzel 1989; Gutkin and Ermentrout 1997. This dynamical approach uses a phase variable φ to model the location along the spike trajectory in the phase plane. The phase evolves according to the stochastic differential equation **

*dφ*=[*β*(1 + cos *φ*)+(1 − cos *φ*)] *dt*+ *σ*(1 + cos *φ*) *dW*t

**for 0 ≤ φ≤ 2π. Here, dWt designates white noise, σ is the noise standard deviation, and βis the bias parameter which models the strength of the input current; a spike occurs every time φ passes through π. In the simulations we used β = σ = 1. Table 2 contains summary statistics for these two simulated spike trains. In each spike train, there are about 6375 spikes, the ISI means are around 3.1, and the ISI histograms are slightly skewed to the right. These two data sets have low CVs, both around 0.33. Figure 2 shows a smoothed histogram for one of the data sets, using a 0.06-ms binwidth. **

**The physiological data are maintained discharges of retinal ganglion cells of goldﬁsh. These data were provided by Dr. Michael Levine of the University of Illinois; he has detailed his laboratory methods elsewhere (Levine 1991; Levine et al. 1988), so we give a brief summary here. The retinae were removed from the goldﬁsh, and maintained in a ﬂow of moist oxygen. Recordings of the retinal ganglion cells were made with an extracellular microelectrode. A level detector signaled the nerve impulses, whose time of occurrence was recorded to an accuracy of within 1 ms. Continuous blocks of data were recorded for at least 20 s in the absence of any change in the illumination of the retina. When the illumination was changed, enough time was allowed so that the behavior appeared stationary. For this paper, we have used data from ten different cells, labeled 68, 72, 75, 78, 80, 82, 83, 87, 90, and 98; Table 3 contains summary statistics for each cell. The recording time for each cell was 30 s. The number of spikes ranges from 703 to 1028. The ISI means range from 29.1 to 42.7 ms. The ISI histograms are all skewed to the right, and the CVs all exceed 1. Figures 3 and 4 show representative smoothed histograms for the data from cells 75 and 78, respectively. **

**We compare the following two models for the ISI histograms: the lognormal and generalized inverse Gaussian. For y ≥ 0, the lognormal has density **

**1 (log y− α1)^{2 }**

*f*(*y**|**α*)= _{y }√ _{2πα}_{2 }exp − _{2α}_{2 }(1)

**where −∞ <α1 <∞ and α2 >0. The generalized inverse Gaussian has density **

*h*(*y**|**ψ,χ,λ*)= ^{(ψ/χ}�^{)}^{λ/2 }y ^{λ−1 }exp ^{� }− ^{1}(*ψy*+ ^{χ }) ^{� }, (2) 2*K**λ*( *ψχ*)^{2 }^{y }

**where Kλ is the modiﬁed Bessel function of the third kind with index λ (Abramowitz and Stegun 1965), **

^{� }∞ ��

*K**λ*(*u*)= ^{1 }x ^{λ−1 }exp − ^{u }(*x*+ ^{1}*) dx *

**22 x**

**0 **

**Let θ =( ψ,χ,λ), and let Θ be the parameter space. Then for −∞ <λ<**

** **

*∞*, Θ = *Θ**λ*, where

**λ **

*{*(*ψ,χ*): *ψ>*0*,χ*≥ 0} if *λ>*0 *Θ*λ = *{*(*ψ,χ*): *ψ>*0*,χ>*0} if *λ*=0

*{*(*ψ,χ*): ψ ≥ 0*,χ>*0} if *λ<*0.

**It is convenient to reparametrize this model thus: ω = ψχ and η = **

*χ/ψ; η is a scale parameter and ω is a concentration parameter when ω> 0. When ω = 0, these interpretations are no longer valid. In that case, it is more appropriate to regard (2) as a gamma or reciprocal gamma density; the norming constants are obtained from the asymptotic formula for λ>0 and ω ↓ 0: *

*Γ*(*λ*)2^{λ−1 }

*K**λ*(*ω*) ∼

*ωλ *

**and the relation Kλ(ω)= K−λ(ω). Below, we express the test statistics in terms of the estimates of β =(ω,η,λ) rather than θ =(ψ,χ,λ). **

**The lognormal model has been used before to ﬁt spike train data from the cat cerebral cortex and other preparations (Tuckwell 1988). Also, Levine (1991) compared the lognormal with the hyperbolic normal, hyperbolic gamma, and inverse Gaussian [which has λ = − ^{1}_{2 }in (2) above] using ad hoc curve ﬁtting methods. He argued that the lognormal was the preferred model among them. On the other hand, the generalized inverse Gaussian family has only sporadically been applied in the neuroscience context. Jorgensen (1982) used it to study a single spike train, and L´ansk´**

**y et al. (1990) suggested its use, but a systematic study of this family in this context has not been done. **

**We propose the generalized inverse Gaussian model for several reasons. First, diffusions are commonly used to model the underlying processes of neurons which lead to their ﬁring. For λ ≤ 0, Barndorff-Nielsen et al. (1978) have shown that the generalized inverse Gaussian distribution arises as the ﬁrst passage time of a diffusion (with certain inﬁnitesimal drift and variance parameters) to a constant boundary. Thus, the estimates of the parameters will inform us about possible mechanisms leading to the spike train. Next, it contains as special cases models that have been proposed for ISI distributions, such as the exponential, gamma, reciprocal gamma, and inverse Gaussian; thus, it allows for a uniﬁed approach to the analysis of ISI histograms. Finally, the generalized inverse Gaussian has the following statistical properties (Jorgensen 1982): It is an exponential family, so it is an analytically convenient and ﬂexible family for model selection. All its members are unimodal with a positive mode unless χ=0 and0 <λ≤ 1, in which case the density is inﬁnite at zero. The range of CVs for this family spans all positive numbers; the correspondence is not one-to-one, as many values of (ψ,χ,λ) give rise to the same CV. Based on extensive numerical work, it is conjectured, but not yet proven, that all members of this family have positive skewness and kurtosis. **

**In this section, we describe the diffusions that give rise to the generalized inverse Gaussian family, and then discuss their potential neurophysiological signiﬁcance. Suppose that θ =( ψ,χ,λ) ∈ Θ, and λ ≤ 0 (so that χ> 0). For x> 0, let σ^{2}(x) be a strictly positive differentiable function, and suppose that **

**x **

*du *

*γ*(*x*)=

*σ*(*u*)

**0 **

**is ﬁnite for all x>0, and that γ(x) →∞ as x→∞. Also, let **

291 | |||
---|---|---|---|

Table 3. Summary statistics for the ten goldﬁsh retinal ganglion cells Cell: 68 72 75 78 80 82 83 87 Spikes: 814 811 823 975 984 1019 1028 1020 Mean: 36.9 36.9 36.4 30.8 30.4 29.4 29.1 29.4 Median: 15.0 18.0 16.0 10.0 10.0 12.0 12.0 13.0 SD: 39.6 37.6 40.7 38.8 37.1 33.1 33.5 31.3 CV: 1.1 1.0 1.1 1.3 1.2 1.1 1.2 1.1 | 90 998 30.0 13.0 32.9 1.1 | 98 703 42.7 20.0 44.0 1.0 |

**⎧� �� 2 λ− 1 ^{� }Kλ−1(γ(x) ψ)1 d **

*−**σ*(*x*)+ ψ � + *σ*^{2}(*x*)

**⎨ 2 γ(x) _{K}_{λ}_{(γ(x) ψ)}4 dx **

*µ*(*x*)= for *ψ>*0*,λ*≤ 0*,*2*λ*+1 1 d

*σ*(*x*)+ *σ*^{2}(*x*)

**⎩ 2 γ(x)4 dx**

for ψ =0,λ<0

**The result of Barndorff-Nielsen et al. (1978) is the following: **

**Theorem. **Let *{**X*(*t*): t ≥ 0} be a time-homogeneous diffusion process with state space [0*,**∞*), inﬁnitesimal variance *σ*^{2}(*x*), and inﬁnitesimal mean *µ*(*x*) given above. Let *x*0 > 0 satisfy *γ*^{2}(*x*0)= *χ*.If λ ≤ 0, then the ﬁrst passage time T = inf*{**t*> 0: *X*(*t*)=0} has the generalized inverse Gaussian distribution with parameter (*ψ,χ,λ*).

**The ordinary inverse Gaussian was ﬁrst used by Gerstein and Mandelbrot (1964); its parameters ( ψ,χ,− ^{1 }) correspond to σ^{2}(x) ≡ σ^{2},**

**� 2 **

**√ **

*µ*(*x*) *≡−**σψ*, and *x*0= *−**σχ*. This elementary example emphasizes the fact that many choices of *σ*^{2}(*x*) yield the same ﬁrst passage time distribution. This non-uniqueness allows for a wide range of diffusions to be covered by this theorem. For instance, Barndorff-Nielsen et al. (1978) showed that the reciprocal gamma arises from two diffusions with distinct inﬁnitesimal parameters: one with constant *σ*^{2}(*x*) and *µ*(*x*) proportional to *x*^{−1}, and another with σ^{2}(*x*) a power of x and *µ*(*x*) ≡ 0.

**Given any diffusion {X(t): t> 0} with inﬁnitesimal parameters µ(x) and σ^{2}(x), we now describe a general method of constructing other diffusions leading to same ﬁrst passage time distribution. For any function g(x) that is strictly monotone and differentiable for x> 0, let X˜(t)= g(X(t)) be a transformed process. Its state space is the range of g, and its inﬁnitesimal mean and variance are **

*µ*˜(*y*)= ^{1 }*σ*^{2}(*x*)g ^{��}(*x*)+ *µ*(*x*)g ^{�}(*x*)

**2 **

**and **

*σ*˜^{2}(*y*)= *σ*^{2}(*x*)[g ^{�}(*x*)]^{2 }

**respectively, where y = g(x) (Karlin and Taylor 1981). Since g is strictly monotone, **

**T = inf {t>0: X(t)=0} = inf{t>0: X˜(t)= g(0)} **

**when X(0) = x0 and X˜(0) = g(x0), the same generalized inverse Gaussian distribution arises from two distinct inﬁnitesimal mean and variance functions. **

**For example, the values λ = − ^{1}_{2 }and ψ = 0 can lead to σ^{2}(x) ≡ 1 and µ(x)=0for x> 0. Applying the transformation g(x)=e^{−x }, we have σ˜^{2}(y)= y^{2 }and µ˜(y)= y/2 for the transformed state space 0 <y ≤ 1. Thus, the ordinary inverse Gaussian distribution is applicable to this diffusion with more complicated inﬁnitesimal parameters. Several models of neural ﬁring that include reversal potentials or shunting inhibition lead to a linear inﬁnitesimal mean, and quadratic inﬁnitesimal variance on a bounded state space (see, for example, Hanson and Tuckwell 1983; L´y et al. 1990). Thus, a subclass of the generalized inverse Gaussian**

**ansk´family arises out of diffusions that have been independently derived from models for neural ﬁring. Of course, some elementary transformations are needed to depict more realistically the state space, which for the Hanson and Tuckwell model is {y : VI <y ≤ VE}, where VI is the inhibitory reversal potential, and VE is the excitatory reversal potential (which exceeds the ﬁring threshold). **

**In many cases, there is not enough information to assume that the candidate families contain the true model that is generating the data. This is certainly the case for the problem at hand, so we leave the true model for the ISI histograms unspeciﬁed; using the data, we then seek to identify the model which is closer to the unknown true model. The approach here is based on the results of Foutz and Srivastava (1977), Nishii (1988) and White (1982). We now describe the approach in general, and then specialize to the case of the lognormal versus generalized inverse Gaussian families at the end of this subsection. **

**To assess the closeness of two densities f and p on the real line, we use the Kullback-Leibler information criterion (or divergence) **

**� ∞ **

*p*__(__*x*__)__

*I*[*p*: *f*]= *p*(*x*) log *d*x

*f*(*x*)

**0 **

**It is known that I[p : f] ≥ 0, with equality only when p = f (Kullback 1968; White 1994). For example, if p(x)= µe^{−µx }and f(x)= λe^{−λx }for x> 0 and zero otherwise, then I[p : f]= λ/µ − 1 − log(λ/µ), which is 0.09 when the ratio of the means of these exponential densities is λ/µ =1.5. Next, if p and f are normal densities with means µ and ν, respectively, and both variances equal to 1, then I[p : f]=(µ− ν)^{2}/2, which is 0.125 when |µ− ν| =0.5. **

**The Kullback-Leibler information criterion is not a metric, but it does have an advantage over metrics for the problem of comparing two approximations f1 and f2 to an unknown density p. Speciﬁcally, **

**� ∞ **

*f*__2____(__*x*__)__

*I*[*p*: *f*1] − *I*[*p*: *f*2]= *p*(*x*) log

*f*1(*x*)

**0 **

*f*__2____(__*Y*__) __

*= Ep log (3)*

*f*1(*Y*)

**where Ep refers to expectation with respect to the density p, which generates the random variable Y. In practice, given independent and identically distributed data (Y1,...,Yn) from p, an estimate of this quantity is **

**n1 ^{� }f2(Yi)**

**log **

*nf*1(*Y**i*) 1

**If this estimate is positive, we conclude that f2 is closer to p than f1. Next, let **

**F = {f(y|α): α∈ A} **

**and **

**G = {g(y|β): β ∈ B} **

**be two families of densities under consideration, and let p(y) denote the true unknown density that is generating the data. Let α∗ be the quasi-true parameter in A associated with the true density p(y): f(y|α∗)isthe member of F that minimizes I[p : f(·|α)], over α ∈ A, and deﬁne β∗ ∈ G similarly. Thus, **

**� ∞ **

*p*__(__*y*__)__

*I*F = *I*F [*p*: *α**∗*]= *p*(*y*) log *dy*

*f*(*y**|**α**∗*)

**0 **

**characterizes the proximity of F to p, and IG = IG [p : β∗] (deﬁned similarly) characterizes the proximity of G to p. Thus, the difference **

**� ∞ **

*I*G − *I*F = *p*(*y*) log ^{f(y|α}^{∗}^{) }*dy *

*g*(*y**|**β**∗*)

**0 **

**292 **

**allows for a comparison of F and G , using estimates of α∗ and β∗, which we describe next. **

**Given a random sample Y =( Y1,... ,Yn) generated from p(y), let the quasi-log-likelihood function of α under F be **

*n*

*L**n*(*α*| Y )= log *f*(*Y**i*| *α*)

*i*=1

**and a quasi-maximum likelihood estimator (QMLE) ˆ αn be a value which maximizes Ln(α| Y ) over α ∈ A; for G , deﬁne Ln(β| Y ) and βˆn similarly. **

**A natural estimate of IG − IF is **

**n1 ^{� }f(Yi| αˆn)**

*Tn = log *

*ng*(*Y**i*| *β*ˆ*n*)*i*=1

**which is asymptotically normal: **

**√ d n[Tn − (IG − IF )] −→ N(0,σ^{2}) (4) **

**��2 where the variance σ ^{2}= Ep log f(Y | α∗)/g(Y | β∗) is estimated by **

*n*��_{2 }_{σ}_{2}1 ^{� }*f*(*Y**i*| *α*ˆ*n*)

**ˆn = log (5) **

*ng*(*Y**i*| *β*ˆ*n*)*i*=1

**The regularity conditions that guarantee the existence of the quasi-true parameters α∗ and β∗, the consistency of the QMLEs, the asymptotic normality of Tn, and the consistency of ˆn are described in the papers **

*σ*2 by Foutz and Srivastava (1977), White (1982) and Nishii (1988). These conditions are satisﬁed for the comparison below; we omit the details of that veriﬁcation. The decision about which of the families, F or G , is closer to the unknown true model *p*(*y*), now proceeds in the following way. We ﬁrst pose the null hypothesis *H*0: *I*G = against the alternative hypothesis

^{I}F √

*H*1: *I*G =*/I*F . We then compute the statistic *n**T*n and its estimated standard deviation ˆ*σ**n*. Given a conﬁdence coefﬁcient *l*, we construct the approximate conﬁdence interval for (*I*G − *I*F ),

*σ*ˆn *σ*ˆ*n*

*I**n,*l = *T*n − *c*l √ *,T*n + *c*l √

*nn *

**where cl is the upper 100(1 − l/2)-percent point of the standard normal distribution (if l =0.95, cl =1.96). If 0 ∈ In,l, then we cannot reject the null hypothesis that the two models are equally close to the true model; otherwise we reject the null hypothesis. Furthermore, if In,l ⊆ (0, ∞ ), we conclude that F is closer to the true model; if In,l ⊆ (−∞ , 0), we conclude that G is closer to the true model. **

**We now turn to the speciﬁcs for the lognormal and generalized inverse Gaussian families. Let F and G refer to the lognormal and generalized inverse Gaussian family, respectively. The QMLEs ˆ αn =(ˆα1n,αˆ2n)ofthe lognormal parameters are **

*n*

**1 **

*α*ˆ1n = log *Y*i

**n **

**1 **

**and **

*n*

**1 **

*α*ˆ2n = [log *Y*i − *α*ˆ1*n*]^{2 }

**n **

**1 **

**respectively. The QMLEs βˆn =(ˆηn,ωˆn,λˆn) of the generalized inverse Gaussian parameters are not analytically tractable, so we used IMSL (Visual Numerics 1994) to evaluate them numerically, by solving the likelihood equations described in Jorgensen (1982). Let **

*n n*

**1 11 **

*Y*¯= *Y*i and *Y*˜=

*n nYi *

**1 1 **

**Then, Table 4. Comparison of the two runs of the integrate-and-ﬁre model and QMLEs for the generalized inverse Gaussian model **

**Run: 1 **

**2 **

*Comparison of the two families*

**√ **

*n**T*n − 0*.*33 − 0*.*45 *σ*ˆn 0.14 0.11

*QMLEs for the generalized inverse Gaussian model *

*η*ˆn 46.66 51.18 *ω*ˆn 2.44 2.78 *λ*ˆn − 2*.*26 − 2*.*87

*T*n = *λ*ˆn log ˆ*η*n + log[2*K*_{λ}ˆ_{n }(ˆ*ω**n*)] + ^{ωˆ}^{n }[ˆ*η*n ^{−1}*Y*¯+ˆ*η**n**Y*˜] �� � ^{2}_{1 }− *λ*ˆ*n**α*ˆ1n − log 2*πα*ˆ2n − (6)

**2 **

**and **

**� _{�}�2 σ2ˆ**

**ˆn = λn log ˆηn + log(2K_{λ}ˆ_{n }(ˆωn)) − log 2παˆ2n (7) **

**+***λ*ˆn log ˆ*η*n + log(2*K*_{λ}ˆ_{n }(ˆ*ω**n*)) − log 2*πα*ˆ2n**+****(ˆ***η*^{−1}*Y*i +ˆ*η**n*Y^{−1})*−*−*λ*ˆn log*Y*i

**n 2 ^{n}^{i }2ˆα2n**

i=1

**Table 4 contains the results of the comparisons for the two runs of the integrate-and-ﬁre model. For both cases the estimates of ( IG − IF ) are approximately 0.01, and the 95% conﬁdence intervals are in the interval (−∞, 0); hence, we conclude that for both runs the generalized inverse Gaussian family is closer than the lognormal family to the true model. Table 4 also contains the QMLEs for the generalized inverse Gaussian model. For both runs, λ^{ˆ}n < 0, so the Theorem above suggests that diffusion processes are useful approximations for the integrate-and-ﬁre model’s subthreshold behavior. For both runs, the estimated value of λ is less than −2, so that the ordinary inverse Gaussian is not a reasonable ﬁt for either run. Figure 1 depicts a smoothed histogram, along with the ﬁtted densities for the ﬁrst run. While both families ﬁt rather well in the right-hand tail (where ISI > 50), the generalized inverse Gaussian is better near the mode. The small value of the estimated difference in the Kullback-Leibler criterion is reﬂected in the rather good ﬁts of both models to the data. The picture for the second run is similar, so we omit it. **

**Similarly, Table 5 contains the results of the comparisons for the two runs of a variant of the Morris-Lecar dynamical model. The estimates of ( IG − IF ) are approximately 0.02; once again, their 95% conﬁdence intervals lead us to conclude that for both runs the generalized inverse Gaussian family is closer than the lognormal family to the true model. Table 5 also contains the QMLEs for the generalized inverse Gaussian model. For both runs, λ^{ˆ}n < 0, so the Theorem above applies once again. However, this time the QMLE of χ is zero, so the special case of the reciprocal gamma is the best ﬁt here; the estimated values of λ are less than −12. Figure 2 depicts a smoothed histogram and **

**Table 5. Comparison of the two runs of the dynamical model and QMLEs for the generalized inverse Gaussian model (reciprocal gamma) **

**Run:1 2 **

*Comparison of the two families*

**√ **

*n**T*n *−*1*.*73 *−*1*.*47 *σ*ˆn 0.14 0.10

*QMLEs for the generalized inverse Gaussian model (reciprocal gamma) *

*χ*ˆn 67.87 71.91 *λ*ˆn *−*11*.*86 *−*12*.*46

the ﬁtted densities for the ﬁrst run of the dynamical model simulation. While neither family reproduces the smoothed histogram, the generalized inverse Gaussian is a better ﬁt at the mode. The extreme tail (ISI > 8) is well ﬁt by both, but the intermediate tail (5 < ISI < 8) is not well ﬁt by either, due to several lesser modes there. The picture for the second run is similar, so we omit it.

Table 6 contains the results of the comparisons for the ten goldﬁsh retinal ganglion cells. Just as for the simulated spike trains, we conclude that for each cell the generalized inverse Gaussian family is closer than the lognormal family to the true model. The estimates of (*I*G − *I*F ) are all between 0.05 and 0.08; once again, their 95% conﬁdence intervals lead us to conclude that for both runs the generalized inverse Gaussian family is closer than the lognormal family to the true model. Table 6 also contains the QMLEs for the generalized inverse Gaussian model. For all cells *λ*^{ˆ}n < 0. In fact, except for cells 72, 78, and 98 the estimated value of λ is near −^{1}_{2}, indicating that the ordinary inverse Gaussian is a reasonable ﬁt for seven of the ten cells. Since the results for the ten cells are otherwise similar, in Figs. 3 and 4 we illustrate cells 75 and 78, respectively, as representatives. The smoothed histograms are considerably noisier than those of the simulations, as evidenced by the many small modes in the 50–150 ms range. The generalized inverse Gaussian is a better ﬁt near the mode, but both models fail to reproduce the tail behavior, since the two families of densities are unimodal. It is unclear at present whether the smaller modes in the tail are the result of noise or indicate some oscillatory activity.

**Fig. 1. Simulated integrate-and-ﬁre spike train **

**The use of diffusion processes to model a neuron’s membrane potential has a long history. Interest in such models is still quite strong, as evidenced by recent work that attempts to understand the mechanisms that lead to highly irregular ﬁring patterns (Softky and Koch 1993). The most thoroughly studied diffusions in this context are the Wiener and Ornstein-Uhlenbeck processes; the latter arises as the continuous limit of Stein’s integrate-and-ﬁre model. The Wiener process yields the inverse Gaussian distribution for the ISI, which is analytically tractable, so that standard inferential procedures (parameter estimation, conﬁdence intervals, goodness of ﬁt tests) have been well developed. That is not yet the case for the Ornstein-Uhlenbeck or other models. This paucity of statistical methods for biologically more realistic models has led to the use of rather ad hoc procedures, such as the use of moments (plots of coefﬁcient of variation versus the mean), or curve ﬁtting using generic models for positive random variables, such as the lognormal. Also, Levine (1991) has suggested that it is preferable to use such procedures because it is hard to estimate separately the various parameters, such as the threshold, leakage time constant, and input rates, of biologically more realistic models. Recent work (Inoue et al. 1995; Iyengar 1997; Iyengar and Liao 1997) has alleviated some of these difﬁculties, but much still remains to be done. **

**In this paper, we have applied a formal inferential procedure for comparing models to assess the adequacy of the ﬁt of the generalized inverse Gaussian family to ISI histograms. For spike trains from simulations and from retinal ganglion cells of goldﬁsh, it is closer, in terms of the Kullback-Leibler information criterion, to the true ISI density than is the log**

**Table 6. Comparison of the ten goldﬁsh retinal ganglion cells and QMLEs for the generalized inverse Gaussian model **

Cell 68 72 75 Comparison of the two families√ nTn −1.61 −1.45 −1.63 ˆσn 0.11 0.09 0.12 | 78 −2.57 0.18 | 80 −2.26 0.12 | 82 −1.90 0.14 | 83 −1.96 0.14 | 87 −1.65 0.13 | 90 −1.89 0.12 | 98 −1.46 0.09 | |||
---|---|---|---|---|---|---|---|---|---|---|

Cell | 68 | 72 | 75 | 78 | 80 | 82 | 83 | 87 | 90 | 98 |

*QMLEs for the generalized inverse Gaussian model *

*η*ˆn 31.00 25.08 31.61 46.07 27.13 32.11 30.68 25.85 27.86 28.10 *ω*ˆn 0.47 0.51 0.45 0.33 0.42 0.47 0.44 0.51 0.46 0.47 *λ*ˆn *−*0*.*36 *−*0*.*17 *−*0*.*39 *−*0*.*76 *−*0*.*40 *−*0*.*57 *−*0*.*54 *−*0*.*39 *−*0*.*44 *−*0*.*16

**Fig. 2. Simulated spike train from the dynamical model **

**Fig. 3. Spike train for goldﬁsh retinal ganglion cell 75 **

**Fig. 4. Spike train for goldﬁsh retinal ganglion cell 78 **

normal. The ﬁt of this model to the simulated data is better than that for the actual spike trains, largely because there are several modes in the right-hand tail of the latters’ ISI histograms. The source of these modes is not clear. They could be the result of further noise that is not accounted for, or of some oscillatory (periodic) activity. The theoretical justiﬁcation given in Sect. 2.3 points to one possible resolution of this difﬁculty, namely, the use of mixtures of these densities to account for the other modes. In particular, given the inﬁnitesimal variance *σ*^{2}(*x*) the parameter *χ*, determines the initial point, which might be the resting potential, or some other intrinsic feature of the neuron; the other two paramters, *λ*and *ψ*determine the inﬁnitesimal mean, which largely reﬂects the inputs to the neuron. Thus, a mixture over *χ*could model variations in the reset voltage. We are presently developing ﬁtting procedures for such mixtures.

The theoretical justiﬁcation above also motivates a more detailed study of the underlying diffusions themselves, with an eye towards more detailed modeling of physiological quantities. For instance, we are studying the problem of ﬁtting parameters of the transformed diffusion above to the Hanson and Tuckwell (1983) model. Also, while the case λ = *−*^{1}_{2 }ﬁts many of the ganglion cells, other cases yield estimates of λ that are much less than that. Approximations to these diffusions for the limiting case *λ**→−∞*may yield simpler subfamilies. Finally, there have been numerous reports that the exponential and gamma densities provide good ﬁts for other ISI histograms, and that (for certain limiting parameter values) these simple densities approximate the Ornstein-Uhlenbeck ﬁrst passage time density. Thus, members of the generalized inverse Gaussian family with *λ>*0 may provide good approximations to other diffusions not covered by the result of Barndorff-Nielsen et al. (1978).

*Acknowledgements*. The authors thank Dr. Michael Levine of the University of Illinois at Chicago for the spike train data, and Boris Gutkin of the University of Pittsburgh for the simulations of the variant of the Morris-Lecar model. We also thank an anonymous referee for comments on an earlier draft of this paper. S.I. and Q.L. were supported by Ofﬁce of Naval Research grants N00014-89-J-1496 and N00014-93-1-1064. Q.L. was also supported by the National Institute of Health grants AI28076 and AI07358.

**Abramowitz M, Stegun IA (1965) Handbook of mathematical functions. **

**New York, Dover Barndorff-Nielsen O, Blaesild P, Halgreen C (1978) First hitting time mod**

**els for the generalized inverse Gaussian distribution. Stoch Processes **

**Appl 7:49–54 **

**Ermentrout GB (1996) Type I membranes, phase resetting curves, and synchrony. Neural Comput 8: 979–1001 **

**Ermentrout GB, Rinzel J (1989) Analysis of neural excitability and oscillations. In: Koch C, Segev I (eds) Methods in neuronal modeling. MIT Press, Cambridge, Mass, pp 135–169 **

**Foutz RV, Srivastava RC (1977) The performance of the likelihood ratio test when the model is incorrect. Ann Statist 5: 1183–1194 **

**Gerstein G, Mandelbrot B (1964) Random walk models for the spike activity of a single Neuron. Biophys J 4: 41–68 **

**Gutkin B, Ermentrout GB (1997) Type I membrane excitability as a mechanism for cortical spike train statistics. Technical report, Mathematics Department, University of Pittsburgh **

**Hanson FB, Tuckwell HC (1983) Diffusion approximations to neuronal activity using synaptic reversal potentials. J Theor Neurobiol 2: 127– 153 **

**Inoue J, Sato S, Ricciardi L (1995) On the parameter estimation for diffusion models of single neurons. Biol Cybern 73: 209–221. **

**Iyengar S (1997) Maximum likelihood estimation for the Ornstein-Uhlenbeck process. Technical report, Statistics Department, University of Pittsburgh **

**Iyengar S, Liao Q (1997) Model selection for neural spike train data. Technical report, Statistics Department, University of Pittsburgh **

**Jorgensen B (1982) Statistical properties of the generalized inverse Gaussian distribution. Springer, Berlin Heidelberg New York **

**Karlin S, Taylor HM (1981) A second course in stochastic processes. Academic Press, New York **

**Kullback S (1968) Information theory and statistics. Wiley, New York **

**L´ansk´y P, Smith CE, Ricciardi L (1990) One-dimensional stochastic diffusion models of neuronal activity and related ﬁrst passage time problems. Trends Biol Cybern 1: 153–162 **

**Levine MW (1991) The distribution of intervals between neural impulses in the maintained discharges of retinal ganglion cells. Biol Cybern 65: 459–467 **

**Levine MW, Saleh EJ, Yarnold P (1988) Statistical properties of the maintained discharge of chemically isolated ganglion cells in goldﬁsh retina. Vis Neurosci 1: 31–46 **

**Morris C, Lecar H (1981) Voltage oscillations in the barnacle giant muscle ﬁber. Biophys J 35: 193–213 **

**Nishii R (1988) Maximum likelihood principle and model selection when the true model is unspeciﬁed. In: Rao CR, Rao MM (eds) Multivariate statistics and probability. Academic Press, New York, pp 392–403 **

**Softky WR, Koch C (1993) The highly irregular ﬁring of cortical cells is inconsistent with temporal integration of random EPSPs. J Neurosci 13(1): 334–349 **

**Visual Numerics (1994) IMSL Math/Library. Visual Numerics, Houston, Texas **

**Tuckwell H (1988) Introduction to theoretical neurobiology, vols 1 and 2. Cambridge University Press, New York **

**White H (1982) Maximum likelihood estimation of misspeciﬁed models. Econometrica 50: 1–26 **

**White H (1994) Estimation, inference, and speciﬁcation analysis. Cambridge University Press, New York **