The procedures of all experiments were approved in advance by the local Swedish Animal Research Ethics Committee (permits M32-09 and M05-12). Adult cats were prepared as previously described (Jorntell and Ekerot, 2002, 2006). Briefly, following an initial anesthesia with propofol (Diprivan ® Zeneca Ltd, Macclesfield Cheshire, UK), the animals were decerebrated at the intercollicular level and the anesthesia was discontinued. The animals were artificially ventilated and the end-expiratory CO₂, blood pressure and rectal temperature were continuously monitored and maintained within physiological limits. Mounting in a stereotaxic frame, drainage of cerebrospinal fluid, pneumothorax and clamping the spinal processes of a few cervical and lumbar vertebral bodies served to increase the mechanical stability of the preparation. To verify that the animal was decerebrated, we made EEG recordings using a silver ball electrode placed on the surface of the superior parietal cortex. Our EEG recordings were characterized by a background of periodic 1–4 Hz oscillatory activity, periodically interrupted by large-amplitude 7–14 Hz spindle oscillations lasting for 0.5 s or more. These forms of EEG activities are normally associated with deep stages of sleep (Niedermayer and Lopes Da Silva, 1993). The pattern of EEG activity and the blood pressure remained stable, also on noxious stimulation, throughout experiments (see also Jorntell and Ekerot, 2006).

Spinal neurons (six separate animals) and cerebellar neurons (eight separate animals) were recorded in different experiments. For spinal recordings, a laminectomy was made to expose spinal segments C5-T1. For cerebellar recordings, the bony cerebellar tentorium was removed to expose the caudal part of the anterior lobe of the intermediate cerebellum, in order to access the forelimb region of the cerebellar C3 zone, where all cerebellar recordings were made. In both cases, a pool of cotton-in-agar was built around the recording regions. The dura was removed and the recording region was covered with paraffin oil to prevent drying of the tissue. In the case of spinal recordings, we also made small openings in the pia in order to facilitate penetration of the patch clamp pipettes. Also for spinal recordings only, the recording region was covered with agarose attached to the bone of the vertebrae in order to dampen tissue movements.

Patch clamp pipettes were pulled to 4–14 MOhm and contained a potassium-gluconate based solution, as previously described (Jorntell and Ekerot, 2006). The spinal neurons were recorded from the lower cervical segments (C6–C8) at depths of 1.8–2.8 mm, corresponding to laminae IV–VII. Since they were located above the level of the motor nuclei that contains the alpha-motorneurons, they were labeled spinal interneurons. The cerebellar cortical neurons were recorded from the superficial parts of the forelimb region of the C3 zone. Purkinje cells were identified by the presence of complex spikes. Molecular layer interneurons (MLints) were recorded from the superficial part of the cortex, visible from the surface, above the first Purkinje cell layer. Spiking neurons in the molecular layer were classified as MLints. Golgi cells were recorded from the granule layer located below the first Purkinje cell layer from the surface. They were classified as Golgi cells on basis of their regular spontaneous firing (Van Dijck et al., 2013).

For recordings of series of interspike intervals, we used a HEKA EPC 800 patch clamp amplifier set to current clamp. As the patch clamp pipettes approached a neuron we released positive pressure and applied gentle suction in order to establish a gigaohm seal before break-in. For all neurons, seal resistance varied between 0.8 and 5 GOhm. Suction and sometimes current command steps of one 1 nA were applied to obtain intracellular access. Access resistance was compensated for off-line and was usually between 20 and 80 MOhm. To change the firing frequency of a neuron we used either steady current injections lasting up to several minutes or current step commands lasting 500–4000 ms. For current step commands, when the neuron's stationary state (see below) bridged the time gap between two successive current step commands, the spike firing of the consecutive current step commands was treated as a single stationary state, omitting the interspike interval bridging the steps. The signal was converted to a digital signal using the analog-to-digital converter Power 1401 mkII from Cambridge Electronic Design (CED, Cambridge, UK). The neural responses were sampled at 100 KHz and recorded continuously with the software Spike 2 from CED. Spike2 was also used to identify the spikes and to obtain data of the spike time series. Power spectrograms of the spike firing frequencies were calculated using the multitaper method as implemented in the Chronux toolbox for Matlab (http://chronux.org/).

All raw spike trains used to construct the models were separated into shorter spike trains that were 50 spikes long, with a few exceptions where the original spike train was shorter than 50 spikes but longer than 40. Only trains that were found to have ISI distributions that were stationary and log-normally distributed were used further (see statistics section for details). All log-normal distributions were constructed as maximum likelihood estimates of the distributions of ISIs using the LOGNFIT method in MATLAB (MathWorks). The normalized input of each stationary region was then estimated by x = log(i.s.d). This value together with the firing frequency during the state was used to fit the frequency curve Equation (6b) by iterative least squares estimation using the MATLAB method NLINFIT. The normalized input was then estimated once again using the inverse of Equation (6b), x = log [exp (E⁻¹/c_x) − 1] − Δ_x. The new estimate was used to evaluate the whether the means and standard deviation the models provided could describe the ISI distribution of the region using log-normal, gamma or inverse Gaussian probability density functions (PDFs). The fraction of regions that could not be rejected was seen as the model accuracy.

One reason for using the log-normal distribution over other commonly used distributions is the possibility to use tests for normality in order to select which spike-trains should be used to fit the model. We use the Shapiro-Wilk (SW) test (0.05 significance level), that have been found to have high statistical power for very short samples (Shapiro and Wilk, 1965). Long spike-trains can hence be divided into shorter trains that are used for the model.

Transient events that could affect the statistics of the spike-trains were avoided by only using regions that were found to be stationary by the Kwiatkowski, Phillips, Schmidt, and Shin (KPSS) (p < 0.05) test. Regions that were rejected by the KPSS-test were not used further.

Where ISIs have to be compared to a specific distribution (e.g., to evaluate the constructed model), we used the 2-sample Anderson-Darling (AD) test (Scholz and Stephens, 1987) instead of the commonly used Kolmogorov-Smirnov (KS) test, since the KS-test underestimates the test statistics whenever the mean and variance of the distribution is estimated from the empirical distribution and it is known to have poor statistical power (Shapiro et al., 1968; Stephens, 1974). The performance of AD is comparable to that of SW when the mean and variance is unknown (Stephens, 1974).

The local variability (LV) (Shinomoto et al., 2005), the third and the forth L-moment in the form of the L-skewness and L-kurtosis was used to compare the performance of the gamma, lognormal and inverse Gaussian distributions. The asymptotic approximations of the L-moments have been found to be reliable with sample sizes of 50 or more (Hosking, 1990). The confidence intervals of the constructed models were created by drawing a total of 50,000 spike trains consisting of 50 spikes each, where the input to the model was varied uniformly over the relevant input range (X ∈ [0, 9]). Linear regression using the Matlab method REGRESS was used to validate the constructed model. The coefficient of determination R² is used as a measure to evaluate the constructed linear regression. From Equation (6a) the line should have the form log(i.s.d.) = p₁ x + p₂ where p₁ = 1 and p₂ = 0. For a valid set of model parameters, the confidence intervals of the parameters of the linear regression should envelope these values. In addition to this a, Durbin-Watson test (p < 0.01) was used to test whether the residuals of the linear regression were uncorrelated using the MATLAB method DWTEST. The SW-test was further used to test whether the residuals were normally distributed (p < 0.01).

Our first aim was to analyze the distribution of the inter-spike intervals (ISIs) across the operative range of firing frequencies that these types of cells have previously been described to display under behavior (Edgley and Lidierth, 1987, 1988; Prut and Fetz, 1999; Pasalar et al., 2006; Takei and Seki, 2010, 2013; Badura et al., 2013). We used constant or episodic current injections to repeatedly modulate the neuron's firing frequency across these different levels, illustrated in Figure 1A for a spinal interneuron. During current injections we obtained spike firing levels during which the series of ISIs were stationary (KPSS-test, p < 0.05). The series of ISIs obtained at each firing level were binned in histograms (Figures 1B,C) to investigate the shape of the distribution of the intervals, which all had the skewed shape seen in the figures. Notably, the crossings of specific arbitrary threshold levels described by the subthreshold spontaneous synaptic activity had a different interval distribution than that of the spikes (Figure 1D), a discrepancy which could be explained by direct observations that the apparent spike threshold varied also at short time intervals (Figure 1E). In particular the long tail of distributions of the threshold crossing events in Figure 1D, which was absent in the ISIs, strongly suggests that the statistics of the spike generation was determined by other factors in addition to the synaptic inputs and the absolute value of the membrane potential. However, as the actual occurrence of a spike would inevitably affect the dynamics of the membrane potential (for example by absolute and relative refractoriness) in a way that we could not completely account for, we did not carry out the analysis described in Figure 1D in a systematic fashion.

The log-normal, the gamma and the inverse Gaussian distributions are typically used to describe skewed shapes of the type found for the ISIs. In this paper we focus on the log-normal distribution since it comes with the advantage that it is tied to an underlying normal distribution, which makes it possible to use tests for normality with high statistical power even with small sample sizes to validate the assumption that the ISI distribution can be described by a log-normal distribution (Shapiro et al., 1968). This is important since also temporary stationary states containing a limited number of spikes can then be used. Furthermore, the investigated neurons are typically participating in movement control of limbs where the phases of the movements have durations of at most a few hundred of milliseconds. Consequently, short spike trains can be useful also in this perspective as the aim of the model is to approximate the behavior of the neurons during natural circumstances. The log-normal distribution has previously been used to fit the spontaneous activity of neocortical neurons in the unanesthetized cat (Burns and Webb, 1976) and the spike output of detailed neuron models (Levine, 1991). As we will show below (Figure 7), the log-normal distribution provide at least as good approximations of the spike time distributions as the gamma and inverse Gaussian distributions.

The shape of the log-normal PDF in Equation (1) is determined by two variables, μ and σ, which corresponds to the mean and standard deviation of the distribution's natural logarithm, respectively.

The log-normal distribution provided a statistically good fit to the series of ISIs obtained at each firing level (Figures 1B,C). We also used the log-normal distribution to describe the spike firing statistics across different firing levels for single neurons (Figure 2). For spinal interneurons (spInt, N = 9, membrane resistance, R_in, = 29 ± 11 MOhm) a total of 564 stationary states with log-normal ISI distributions were tested, corresponding to a mean of 63 ± 54 states per neuron. We also tested Purkinje cells (PCs, N = 3, R_in = 4.2 ± 1.2 MOhm), from which we obtained a total of 85 states (28 ± 21 states per neuron), molecular layer interneurons (MLint, N = 3, R_in = 132, 160, and 177 MOhm) with a total of 34 states (11 ± 8 states per neuron) and Golgi cells (Goc, N = 3, R_in = 30, 35, and 270 MOhm) with a total of 241 states (80 ± 46 states per neuron). Because the mathematical descriptions we apply below required a relatively high number of stationary states under stable conditions and for multiple levels of firing, we restricted our initial analysis to this limited set of neurons. In addition, to provide controls for the specificity of the model described later in the paper, we used the spike firing statistics of an additional 45 spinal interneurons.

The parameters of the log-normal distribution can be translated into the ordinary mean and standard deviation of the distribution through Equation (2): (2)E=exp(μ+12σ2) s.d.=Eexpσ2−1 where E is the mean ISI and s.d. is the standard deviation of the ISIs. By translating the parameters in this way, it is possible to relate them to the type of input frequency-current (f-I) curves that are routinely used to characterize neurons (Mckay and Turner, 2005; Molineux et al., 2005; Zhong et al., 2010). For this reason, we will use the inverse of the mean ISI to obtain the mean frequency E⁻¹ below. We will also use the inverse of the standard deviation (i.s.d.) of the ISIs in order to keep the same unit for both measurements.

Our next aim was to design a model describing the relationship between the input to the neuron and the statistics of its spike firing. The relationship between the input current and the spike firing frequency is frequently reported to be linear (Nowak et al., 2003; Mckay and Turner, 2005; Molineux et al., 2005; Meehan et al., 2010; Zhong et al., 2010). However, as the firing frequency of a neuron approaches zero the linearity inevitably disappears, either by a sharp threshold or by a smoother transition between zero and non-zero firing, also shown theoretically using leaky integrate-and-fire (LIF) neuron models (Fourcaud-Trocme et al., 2003; La Camera et al., 2008). Our results suggest a smoother transition in this non-linear region (Figure 3A), which is in agreement with other in vivo data (Priebe et al., 2004) and in vitro data stimulating the neurons with in vivo like inputs (La Camera et al., 2008), both obtained from neocortical neurons. It is possible that this property is more pronounced in an in vivo setting such as ours with an intense synaptic background activity as compared to in vitro since previous simulation work has suggested that the background level of stochastic synaptic input changes the dynamics of spike firing (Jaeger and Bower, 1999; Gauck and Jaeger, 2000, 2003; Destexhe et al., 2001; Salinas and Sejnowski, 2002; Fellous et al., 2003; Suter and Jaeger, 2004), especially close to the threshold (La Camera et al., 2008). A suitable description of the spike firing statistics should therefore be capable of featuring both a sharp and a smooth threshold in this region. At the high end it should be approximately linear, but also have a transition from the linear region to asymptotically reach 0 as the neuron is sufficiently hyperpolarized (Figure 3B). Here we introduce Equation (3), which satisfy these requirements. In this equation, the parameter c determines the width of the threshold region with c ~ 0 leading to a sharp threshold and higher values lead to a more smooth transition (Figure 3B).

where I denotes the total input current over the membrane (i.e., the sum of synaptic and injected currents), b determines the position of the threshold and a the slope of the linear region. Due to the asymptotic linearity of Equation (3) at high input, it does not capture the high input saturation of the f-I curves obtained from LIF models (Rauch et al., 2003). However, none of the neurons we recorded displayed this behavior within the range of activity studied. It is possible that this behavior would have appeared if the neurons had routinely been driven to even higher levels of activity. However, when tested for a few of our neurons, at these levels we observed a spike break-down (severe reduction in spike amplitude, broadening of the spike) before the appearance of high input saturation.

An interesting feature became evident when we plotted the i.s.d. against the firing frequency. Rather than being linearly related, as would be expected from naïve rescaling of the ISI distribution, this relationship seemed to be best described as an exponential curve (Figure 4A). It is also clear from the two sample neurons in Figure 4A that the relationship was not the same for different neurons. A general exponential relationship between firing frequency and i.s.d. that could approximate both cases in Figure 4A is described by Equation (4):

If the exponential curve appropriately approximates the relationship between the i.s.d. and input, it allows the model to be divided into two consecutive parts. The first part translates the input (e.g., input current) into a normalized dimensionless input x such that i.s.d. = exp x. Since all models must obey this relationship, it can be used as a method to validate the constructed models against the experimental data. By using a logarithmic scale, the exponential curve is transformed into the linear relation shown in Figure 4B, which simplifies visual inspection, but also enables the use of linear regression to validate the constructed model by testing whether the actual data are distributed along this line.

If the normalized input x is used instead of I in Equation (3), after modifying b and c accordingly, the coefficient of variation (CV) of the ISI over a range of input can be computed using Equation (5). Figure 4C illustrates how the CV varies using the different values of the parameters a and c that were also used in Figure 3B.

During periods of high excitation, the CV will approach 0, as the ISIs will be defined by the properties of the neurons refractory period (Softky and Koch, 1993). The behavior during periods of low excitation is not as well defined. The parameter c of Equation (5) determines whether the CV approach 0 (c < 1), diverge (c > 1), or converge toward a fixed value (c = 1). The latter has support from studies of biophysical integrate-and-fire models (Softky and Koch, 1993; Rauch et al., 2003) while the other two options lack support or acceptable interpretations. From this, the parameter c should be set to 1 for the description of the firing statistics to have any meaning as the firing intensity approach 0. The actual asymptotic value of the CV can be calculated by the formula in Table 1.

The complete model can be summarized by the three equations in Equation (6). The first two Equations (6a,b) form a parametric equation between frequency and i.s.d. with the normalized input x as the parameter. In this way the model separated into one part that translates input to normalized input that is determined by the parameters Δ _I and c_I, and one part that describes the relationship between frequency and i.s.d. that is determined by the parameters Δ _x and c_x. The parameters and some related quantities are described in Table 1.

In summary, to find the parameters of the model that best describes the spike generation of a neuron from experimental data, one would need to find the mean and s.d. of the ISIs from several stationary states. These data can then be used to find the optimal parameters for that neuron in Equation (6).

Figure 5A demonstrates the fit of the spinal interneurons recorded to the relationship between input and the i.s.d. predicted by the model. None of the neurons failed the DW-test (p > 0.01) of the residuals and the linear regressions provided good fits (R² = 0.96 ± 0.02; N = 9). In Figure 5B the relationship between input and spike frequency for each recorded spinal interneuron together with the fitted model of the spike generation for each respective neuron are shown. We further tested the accuracy with which the model could predict the distribution of ISIs for each stationary state. Using AD-test, we found that the model could predict the empirical ISI distribution using the lognormal distribution 97 ± 5% (p > 0.01, N = 564) of the states.

Figure 6 illustrates the applicability of the model to the recorded Purkinje cells, molecular layer interneurons and Golgi cells. As for the spinal interneurons, the predicted relationship between input and the i.s.d. of the model provided a good fit also for the Purkinje cells (p > 0.01; DW-test; R² = 0.94 ± 0.01; N = 3). The model accuracy, using AD-test was 94 ± 5% (p > 0.01; N = 85 stationary states). The model provided a good fit also for molecular layer interneurons (p > 0.01; DW-test; R² = 0.99 ± 0.02; N = 3), with a model accuracy of 96 ± 7% (AD-test; p > 0.01; N = 34 states).

For the Golgi cells, the F-test of the residuals from the predicted relationship between input and i.s.d. failed for 2 out of the 3 neurons (p < 0.01; DW-test; R² = 0.86/0.82), whereas it could not be rejected for the remaining neuron (p > 0.01; DW-test; R² = 0.97). The model accuracy for the individual states was still quite high, 91 ± 5% (p > 0.01; AD-test; N = 241).

Finally, the residuals of all the linear regressions could not be rejected as normally distributed (p > 0.01, SW-test) and the 95% confidence intervals of the linear regression parameters all contained the predicted values (p₁ = 1 and p₂ = 0, see Materials and Methods).

We also compared the accuracy and the higher order moments of the lognormal distribution with those obtained using the gamma and inverse Gaussian distributions. As shown in Figure 7, the approximations obtained using the latter two did not provide better predictions of the higher order moments of the empirical distributions. The accuracy (p > 0.01, AD-test) using the different distributions were also comparable. Using the lognormal distribution, 95 ± 5% accuracy was obtained, while the gamma distribution lead to an accuracy of 97 ± 3%, and the inverse Gaussian distribution to an accuracy of 95 ± 5%.

It has previously been shown that it is possible to distinguish the type of neuron from its spontaneous firing using the firing rate and the entropy of the neurons spike train (Van Dijck et al., 2013). Note that the spontaneous firing in this case refers to the spontaneous activity at rest. In Figure 8, the models that were fitted to the neuron types in this study are compared to each other in order to investigate whether it is possible to differentiate the neurons types also outside of their resting state when they are driven to fire at other intensities. The analysis in Figure 8 is limited to the first two moments of the distributions, due to the short length of the spike trains, but the results indicate that there is no apparent relationship between the type of neuron and model behavior. Rather, the differences between neurons of the same type are comparable to the differences between neurons of different types. A consequence of this relationship, and the relationship between the firing frequency and the standard deviation of the ISIs described above, is that the spike firing statistics in the different neurons we explored here are apparently similar (or rather equally dissimilar) across neuron types when their input is modified so that they have similar firing rates.

To rule out that the accuracy of the constructed models were not simply the result of over-fitting of random data or the possibility that any kind of unrelated ISI data would yield similar results, the accuracy of the models fitted to single neurons were compared to models that were fitted to 4, 50, and 100 stationary states randomly selected out of a total of 498 states (only states that could not be rejected as log-normally distributed, SW-test, p > 0.05) obtained from 45 different spinal interneurons with spontaneous activity (>2 Hz).

The proportion of the data points obtained from the randomly selected states that was accurately predicted by the model was on average 67 ± 13% (4 random states: 76 ± 28%, 50 random states: 67 ± 8%, 100 random states: 67 ± 5%; p < 0.01; AD-test), which was significantly lower than for the single neuron models (p < 0.01; student's t-test) for all neuron types. The distribution of R² (4 states: R² = 0.84 ± 0.25; 50 states: R² = 0.88 ± 0.06; 100 states: R² = 0.89 ± 0.02) was, just like the model fit, significantly different compared to the single neuron models (p < 0.01; student's t-test), except for the Golgi cell models.

One requirement to perform the present analysis was to obtain long-lasting intracellular recordings under stable conditions in vivo, for which the whole cell recording technique was our method of choice. But for this approach there are sources of variability for the neuronal activity that could not be avoided, and which could change the relationship between the injected current and the average firing frequency. These sources are variations in background synaptic activity and variations in access resistance between the electrode and the cell. Although access resistance was compensated for off-line, this cannot be done with 100% accuracy and the error can vary over time. The variations in the background activity were outside of our control.

In order to validate the assumed linear relationship between the input current and the normalized input Equation (5c), the model was used to predict the normalized input from the experimental data. For each stationary state, both the s.d. and the frequency of the state were used to predict the normalized input through the inverse of Equations (6a,b). These predictions should be linearly related to the actual input current according to Equation (6c).

Figure 9 illustrates the relationship between injected current and the predicted normalized input for four sample neurons with a high number of stationary and log-normally distributed states with different current injections. The variability of the predicted input for each input current level is due to both the inherent stochasticity of the spike generation, but also the issues with background synaptic noise and variations in access resistance. However, overall there was a relatively good correspondence between the amount of injected current and the estimated normalized input, with the variability lying in a similar range as reported from other in vivo recordings (Nowak et al., 2003).

A model of the firing statistics neurons across its potential range of stationary states can be extended to inhomogeneous situations with input that is not constant. The least complex model that does not incorporate transient behavior would use the log-normal hazard function in Equation (7) to describe the instantaneous firing rate.

where λ is the instantaneous firing rate, f is the log-normal PDF in Equation (1) and Φ is the cumulative density function of the standard normal distribution.

Intracellular recordings of the membrane potential as well as the spike firing during input modulation equivalent to that which would occur under behavior does to our knowledge not exist for any regular spinal interneuron in vivo, but does exist for the spinocerebellar tract cells of the lumbar cord. To demonstrate the applicability of our model in inhomogeneous situations, the parameters adapted to one of our spinal interneurons were used to simulate the firing during slowly modulated input that is seen during fictive locomotion in DSCT neurons (Fedirchuk et al., 2013). The recorded membrane potential, averaged across a number of cycles of locomotion, was used as the input to the modeled neuron (Fedirchuk et al., 2013). Figure 10 illustrates the fit between the output predicted by the model and the recorded instantaneous firing frequency. The actual firing almost completely stayed within the 95% confidence bounds, indicating that under these conditions the inhomogeneous model, generated from a spinal interneuron in the lower cervical segments, accurately predicts the actual firing behavior of DSCT neurons during fictive locomotion (with the exception of one data point of the actual spike firing that fell outside the confidence bounds of the model prediction, see Figure 10). The initial rise of the experimental response indicates a tendency of overshoot following the period of inactivity. Note however, that while the cyclic membrane potential modulations used as input and the instantaneous firing frequency were recorded from the same neuron (Fedirchuk et al., 2013), they were naturally not from the same step cycles. Due to this, some discrepancies between the predicted output by the model and the recorded frequency is expected.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.