The time rescaling theorem states that the interspike intervals (ISIs) of a continuous time point process can be transformed, or rescaled, so that the rescaled process is Poisson with unit rate, e.g. the rescaled ISIs are independent and exponentially distributed. This variable transform takes the form

{

*t*_{i}} is the set of spike times and λ(

*t*|

*H*_{t}) is the conditional intensity function: temporally varying and history dependent spike probability. Although we will henceforth drop the

*H*_{t} in our notation, such conditioning on the previous spiking history is always implied. Intuitively, the ISIs are stretched or shrunk as a function of total spike probability over the ISI interval so that the rescaled ISIs are centered about a mean of 1. Several proofs of this theorem exist, (

Brown, 2001). Here we present a simple proof of the exponential distribution of the rescaled ISIs. A proof of their independence is located in

Appendix A.

The proof proceeds by discretizing time into bins of width Δ, writing down the probability for each discrete ISI, and then taking the continuous time limit, e.g. Δ → *dt*. The discrete time bins are indexed as *k* and the bins within which the spikes occur are denoted as *k*_{i}. Further we define *p*_{k} as the discrete probability mass of a spike in bin *k* and like λ(*t*) it should be taken as conditionally dependent on the previous spiking history.

The probability of the i-th ISI is the probability that there is a spike in bin

*k*_{i} given that the preceding spikes were located in bins

*k*_{1},

*k*_{2},…,

*k*_{i−1}.

where

*L*_{i} is defined such that

*k*_{i−1} +

*L*_{i} =

*k*_{i}. This is simply the product of the probabilities that there are no spikes in bins

*k* = {

*k*_{i−1} + 1,

*k*_{i−1} + 2,…

*k*_{i} − 1} multiplied by the probability that there is a spike in bin

*k* =

*k*_{i}. For simplicity we now drop the

*i* subscripts.

In preparation for taking the small bin size limit we note that when Δ becomes small, so does

*p* e.g.

*p* << 1 for all bins. This implies that 1 −

*p* ≈

*e*^{−p} allowing the above to be rewritten as

Note that the upper limit of the sum has been changed from

*L* − 1 to

*L* with the justification that we are in a regime where all the

*p*’s are very small. We define the lower and upper spiketimes as

*t*_{k} =

*k*Δ and

*t* =

*t*_{k+L} = (

*k* +

*L*)Δ, define λ(

*t*_{k+l}) such that

*p*_{k+l} = λ(

*t*_{k+l})Δ

^{2} and also define the ISI probability

*density* *P*(

*t*) such that

*P*(

*k* +

*L*) =

*P*(

*t*)Δ. After substituting these into the above equation and converting the sum to an integral we obtain

Consulting

equation 1 we note that the integral in the exponent is by definition τ. Further, applying the Fundamental Theorem of Calculus to this integral gives

*d*τ = λ(

*t*)

*dt* ^{3} Changing variables from

*t* to τ we finally obtain

which is now exponentially distributed and completes the proof.

Although the τ

_{i} can be compared to the exponential distribution, it is useful to note that a second variable transform will make the rescaled ISIs uniformly distributed.

General practice is to sort the rescaled ISIs

*z*_{i} into ascending order and plot them along the y axis versus the the uniform grid of values

where

*N* is the number of ISIs and

*i* = 1,…,

*N*. If the rescaled ISIs

*z*_{i} are indeed uniformly distributed then this plot should lie along the 45 degree line. Essentially, the cumulative density function (CDF) of the rescaled ISIs

*z*_{i} is being plotted against the CDF of the uniform distribution (the

*b*_{i}’s). We show an example of such a plot in . Such a plot can be thought of as a visualization of a Kolmogovov Smirnov (KS) test which compares two CDFs and is usually referred to as a “KS plot.” Formally we can state the null hypothesis

*H*_{0} of this test as follows:

*H*_{0} *: Given a model of the conditional intensity function which is statistically adequate, the experimentally recorded ISIs can be rescaled so that they are distributed in the same manner as a Poisson process (exponentially distributed) with unit rate.*

Under the null hypothesis the maximum distance between the two CDFs will, in 95% of cases, be less than

where

*N* is the number of rescaled ISIs (

Brown, 2001;

Johnson and Kotz, 1970). Equivalently, the plotted line of rescaled ISIs will lie within the bounds

in 95% of cases under the null hypothesis. It should be kept in mind that this is

*not* equivalent to saying that the line of rescaled ISIs lying within these bounds implies a 95% chance of the model being correct.

2.1 Temporal Discretization Imposes KS Plot Bias

The time rescaling theorem applies exactly to a point process with instantaneous events (spikes) and infinite temporal precision, e.g. continuous time. As a practical matter when fitting a statistical model one generally discretizes time. For discrete time the integral of

equation 1 is naively replaced by a sum

If

*p*_{k} << 1 ∀

*k*, e.g. in situations where either the bin size is very small Δ → 0 or the firing rate is very low, the time rescaling theorem will apply approximately even if a discrete time model is used. However, it often happens that

*p*_{k} is in fact large. For example 50 Hz spiking sampled at 1 msec implies

*p* ≈ 0.05, and under many conditions the firing rate can be much higher, at least over some subset of the recording, e.g. during bursting. In such cases, the rescaled times τ

_{i} will

*not* be exponentially distributed and the KS plot will exhibit significant biases (divergences from the 45% line) even if the discrete time model for

*p*_{k} is exactly correct. We this demonstrate in where two KS plots generated using the exact same spikes and time varying firing rate are shown, but a temporal discretization was imposed for one of the plots.

These biases originate in two distinct consequences of discretizing a continuous process. First, there is a lower bound upon the smallest possible ISI (one bin) which leads to a lower bound on the smallest possible rescaled time *z*. Second because only a single spike per bin is allowed, using a discrete time model to estimate the firing rate of a continuous time process results in an underestimation of the firing rate. To demonstrate these issues fully we now consider the simple case of a homogeneous Bernoulli process with a constant spike probability *p*_{k} = *p* per bin for which the the CDF of the *z*’s can be calculated analytically and the KS plot determined exactly.

For a discrete time process only a discrete set of ISIs are possible. Specifically {

*n*Δ} where

*n* is an integer greater than zero and Δ is the bin width. In the case of a homogeneous Bernoulli process the rescaled ISIs are τ(

*n*) =

*pn* and

and the discrete probability distribution of interspike interval times (and rescaled times) is

As in

equation 2, this is merely the product of the probability of no spike for

*n* − 1 bins followed by the probability of a spike in the last (

*n*th) bin. The ’B’ subscript indicates the Bernoulli process.

*P*_{B}(

*n*) is not an exponential distribution, as would be expected for a homogeneous Poisson process. It is a geometric distribution, although in the limit of small p, it reduces to an exponential distribution.

^{4} The CDF of this ISI distribution is easily calculated by summing the geometric series and combining terms.

To finally get the CDF of the rescaled ISIs

*z*,

equation 9 is inverted to get

and substituted into

equation 11.

In we use

equation 12 to generate the KS plot for various spike per bin probabilities

*p*. Even at

*p* = 0.04 which would correspond to 40 Hz firing at 1 msec discretization the CDF is highly non-uniform with a step like structure caused by the discrete values which the rescaled ISIs can take. Such “steps” will be smoothed out if a inhomogeneous Bernoulli process is used instead. (see below) There is, however, another more serious divergence from uniformity, namely

*a distinct positive bias at low (close to 0) rescaled ISIs and a distinct negative bias at high (close to 1) rescaled ISIs*. This bias will not disappear if an inhomogeneous Poisson process is used.

The grey band, which is barely visible, is the 95% confidence region of the KS plot assuming 10 minutes of 40 Hz spiking, which translates into 24000 spikes on average. Since the confidence bounds are so close to the 45 degree line, and will be for any spike train with a long recording time and appreciable firing rate, we introduce a new type of plot in that we term a “Differential KS Plot”. This is simply a plot of the difference between the distribution we hypothesize that the CDF of rescaled times should follow (in this case uniform) and the CDF of the “experimentally recorded” rescaled ISIs (in this case the rescaled ISIs of the Bernoulli process). E.g.

It should be emphasized that the Differential KS plot displays the same information as the KS plot, but does so in a different and more visually accessible manner. The confidence bounds (horizontal dashed lines) are now simply given by

where

*N* is again the number of rescaled ISIs. Plotted this way one can clearly see the positive bias at low values of the rescaled ISIs and the negative bias at high values of the rescaled ISIs. We emphasize that since these KS and Differential KS plots are calculated using the exact homogeneous Bernoulli distribution, the biases are not finite sampling effects.

The positive bias at low ISIs is easily understood by noting that the smallest possible rescaled time is not zero but

What about the negative bias at large (z close to 1) rescaled ISIs ? Consider a homogeneous Poisson process with a firing rate λ. Upon discretizing time into bins of width Δ one might naively expect the probability of a spike per bin to be

*p* = λΔ. However it is in fact slightly less than this as we now show. Assume that there is a spike at time

*t* = 0. Then for a homogeneous Poisson process the probability density for the waiting time

*t*_{w} until the next spike is ρ(

*t*_{w}) = λ

*e*^{−λtw}. Integrating, the probability that the next spike lies within any interval

*t* <

*t*_{w} ≤

*t* + Δ can be obtained.

Defining the bin index

*n* such that

*t* = (

*n* − 1)Δ and discretizing we get

where we have defined

*p* = 1 −

*e*^{−λΔ} in the last line. Discretizing time transforms the homogeneous Poisson process into a homogeneous Bernoulli process, but with a per bin probability of a spike

*p* ≠ λΔ. In fact, by expanding the exponential as a Taylor series it can be seen that

The continuous Poisson process still has an expected number of spikes per interval of width Δ of

, but such an interval could have more than one spike in it. In contrast, the discrete Bernoulli process can only have either 0 or 1 spikes per bin. Therefore the per bin “spike probability”

*p* calculated above is not the expected number of spikes of the continuous point process within an interval Δ. It is the expected number of

*first* spikes in an interval Δ, which is of course less than the total number of expected spikes. By discretizing into bins any chance of there being more than one spike in a time window Δ has been eliminated.

The breakdown of the first order expansion of the exponent is the source of the negative KS plot bias at high (z close to 1) rescaled ISIs. It is a fundamental consequence of discretizing a continuous time point process and is closely connected how the conditional intensity function is generally defined, that is, as the small bin size limit of a counting process (see (

Snyder, 1975) cf. Chapter 5). More specifically the conditional intensity function is the probability density of single spike in an infitesimal interval [

*t*,

*t*+Δ). As shown above, this probability density is actually

*p*/Δ = (1−

*e*^{−λΔ})/Δ < λ and the equality only holds in the limit. Thus

*p*/Δ is not a good approximation for λ when the bin size is too large, and this causes the time rescaling theorem to break down.

2.2 Inhomogeneous Bernoulli Processes

The same positive (negative) bias in the KS plot at low (high) rescaled ISIs remains when the spiking process is not homogeneous Bernoulli. We now we define three inhomogeneous spiking models in continuous time and subsequently discretize them. We use these inhomogeneous discrete time models to simulate spikes and then calculate the rescaled ISIs using the exact discrete time model used to generate the spikes in the first place. The goal is to show that even if the exact discrete time generative model is known, the continuous time rescaling theorem can fail for sufficiently coarse discretizations.

The first model is an inhomogeneous Bernoulli process. One second of the inhomogeneous firing probability is shown in . The specific functional form was spline based with knots spaced every 50 msec and the spline basis function coefficients chosen randomly. This model firing probability was repeated 600 times for a total of 10 minutes of simulated time. The second and third models were the homogeneous and inhomogeneous Bernoulli models respectively but with the addition of a spike history dependent renewal process shown in . We used a multiplicative model for the history dependent firing probabilities of the form

where λ

_{0}(

*t*) is the time dependent firing probability independent of spike history effects and λ

_{hist} is the spike history dependent term which is a function of the time since the last spike (

*t*′ =

*t* −

*t*_{ls}). The functional form of the spike history dependent term was a renewal process, specifically

where

*t*′ =

*t* −

*t*_{ls} is in msec. This form was chosen to mimic a brief refractory period and subsequent rebound. For comparison purposes, all three of these models were constructed so that the mean firing rate remained approximately 40 Hz. Thus the inhomogeneous Bernoulli firing probability had a 40 Hz mean. In the spike history dependent cases, the history independent firing probabilities λ

_{0}(

*t*) were adjusted downwards so that when history effects were included the mean firing rate remained approximately 40 Hz. Specifically the history independent firing probability of the homogeneous Bernoulli process was reduced to 29 Hz and a similar reduction was made for the inhomogeneous Bernoulli model.

In we demonstrate the effect on the KS and Differential KS plots when these models are subjected to various temporal discretizations. Specifically we discretized the models at 1, 0.5, and 0.1 msec resolution by averaging λ

_{0}(

*t*) over these bin widths, e.g.

*p*_{k,0} =< λ

_{0}(

*t*) >

_{k}. The spike history dependent term is a function of

*t*′ =

*t* −

*t*_{ls} which was also partitioned into bins. Similar averaging was then employed so that

*p*_{k′,hist} =< λ

_{hist}(

*t*′) >

_{k′}. The full discrete conditional spike probability is then

*p*_{k} =

*p*_{k,0}*p*_{k−kls,hist} where

*k*_{ls} is defined as the most recent bin prior to bin

*k* that has a spike in it. We then simulated 10 minutes worth of spikes for each model and discretization.

^{5} After generating the spikes we then calculated the rescaled times and CDF difference plots according to

show the results for the inhomogeneous Bernoulli model. Comparison with reveals that the main effect of inhomogeneity is to smooth out the “steps”. The positive (negative) biases at low (high) rescaled times remain, and as expected they are smaller for finer temporal discretizations. show the results when the spike history dependent term is added to both the homogeneous and inhomogeneous Bernoulli models. The important point is that the biases are

*worse* for both models when spike history effects are included, even though the models are constructed so that the mean firing rate remains 40 Hz. The reason for this is that the history dependent term is constructed so that the spike train exhibits burst like behavior. Specifically, after a short (2 msec) refractory period there is an increased probability of a spike. This increases the number of short ISIs. It also increases the smallest possible rescaled ISI

*z* because the probability of a spike goes up immediately following a prior spike, and this ends up shifting distributional weight to short ISIs. This is an important point because it implies that in real experimentally recorded spike trains, which may exhibit burst type behavior, the bias in the KS plot will be worse than would be expected by a simple estimate based upon the mean firing rate, as given in

equation 14.

2.3 Unbiased Discrete Time Rescaling Test Using Model Simulation

In a previous section we showed analytically that when discrete time models are used, the rescaled ISIs may not be exponentially distributed even if the model is exactly correct and that this manifests in the KS plot as systematic biases. Our first proposed solution (we present a second in the following section) to the bias problem is to not assume that the rescaled ISIs are exponentially (or uniformly) distributed, but to instead use a procedure similar to bootstrapping. This proceeds by noting that if a candidate model accurately describes the recorded spikes, then the rescaled ISI distribution of the spikes and *the rescaled ISI distribution expected by the fitted model* should be statistically indistinguishable. If instead the model form is inappropriate to describe the spiking data, then the rescaled ISI distribution expected by the candidate model will not match that of the experimentally recorded spikes, because the model does not describe the recorded spikes accurately. Although the expected distribution of rescaled ISIs is implicitly defined by the fitted model, in practice an explicit analytic form for this distribution may be hard to come by. It can however be sampled numerically using the fitted model to generate spikes and rescaling the resulting ISIs as a function of the model used to generate them. ^{6}

Specifically, after a candidate model is proposed and fit to the recorded spike train data (any type of candidate model may be used as long as it provides an estimate of the conditional intensity function λ) we use the model to simulate spikes, rescale the resulting simulated ISIs and then use a two sample Kolmororov Smirnov test to determine if the sample of estimated rescaled ISIs {

*z*_{est}} and the sample of experimentally recorded rescaled ISIs {

*z*_{exp}} are consistent with being drawn from the same underlying distribution (

Press, 2007). Formally, the null hypothesis of the KS test has been changed from that stated in section 2 and adapted to the case of discrete time data using a model based approach.

*H*_{0}(*estimated*) *: Given a model of the conditional intensity function which is statistically adequate, the set of experimentally recorded ISIs can be rescaled so that they are distributed in the same manner as a sample of rescaled ISIs generated by the statistical model itself.*

To determine the number of spikes (or length of time) which must be simulated we use the analytical expression for the confidence bounds of the two sample KS test. These bounds are a function of the sample sizes of the distributions being compared, in our case the size of the empirical distribution

*N*_{exp} dictated by experiment and the size of the simulated distribution

*N*_{sim} which we can chose. The two sample KS test determines the maximum difference between the CDFs of the experimentally recorded rescaled ISIs {

*z*_{exp}} and the set of rescaled ISIs simulated using the model {

*z*_{sim}}. If the maximum difference between the two CDFs is less than a certain value, specifically

then the null hypothesis is confirmed at the α = 0.05 significance level (

Press, 2007). Alternatively a Differential KS plot (as discussed in the prior subsection) will have 95% confidence bounds of

where we have written

*N*_{sim} = γ

*N*_{exp}.

Since

*N*_{exp} is fixed by experiment, the test will be most strict (tightest confidence bounds) when

*N*_{sim} → ∞, or equivalently as γ increases. Formally, increasing

*N*_{sim} increases the power of the KS test and reduces the number of false positives, e.g. false rejections of the Null Hypothesis. Fortunately

*N*_{sim} need not be overly large. Already at γ = 20 the confidence bounds are only a factor of 1.02 wider than they would be in the infinite limit

in which the exact distribution would be known. This implies that a simulated spike train 20 times longer than the original, experimentally recorded, spike train provides a test power close to optimal and is sufficient to approximate the confidence bounds extremely well. In the results section of this paper we use simulated spike trains 100 times the original “experimental” length (γ = 100) which widens the confidence bounds by a factor of only 1.0005.

Specifically then, this technique proceeds as follows:

Procedure for Numerical Correction - Discretize the spike train into bins of width Δ and fit a discrete time statistical model.
- Rescale the experimentally recorded ISIs using equation 8 to obtain the set of rescaled ISIs {
*z*_{exp}}. - Use the statistical model to estimate the rescaled ISI distribution. Simulate γ > 20 spike trains of the same length as the original experimentally recorded spike train, and rescale the simulated ISIs {
*z*_{sim}}. - Construct the CDFs of both
*z*_{exp} and *z*_{sim}, take the difference, and plot it on the interval [0,1] with the confidence bounds .

2.4 Discrete Time Version of Time Rescaling Theorem

We now prove a discrete time version of the time rescaling theorem which corrects for both sources of KS plot bias. Specifically, we demonstrate how to write a new rescaled time ξ which is exponentially distributed at for arbitrary temporal discretization. The proof given here assumes that an underlying continuous time point process λ(*t*|*H*_{t}) is sampled at finite resolution Δ.

Proposition Suppose a continuous time point process λ(

*t*|

*H*_{t}) is sampled at finite resolution so that the observation interval from (0|

*T*] is partitioned into bins of width Δ. Denote the bin in which the i’th spike is located as

*k*_{i} and that of the next spike as bin

*k*_{i+1} =

*k*_{i} +

*L*_{i}. Let

*p*_{ki+l} =

*p*(

*k*_{i} +

*l*|

*H*_{ki+l}) be the

*discrete time* conditional spike probabilities evaluated in bins

*k*_{i+l} for

*l* = 1…

*L*_{i}. Define the random variable

where

and δ

_{i} [0, Δ] is a random variable determined by first drawing a uniform random variable

*r*_{i} [0, 1] and then calculating

Then ξ_{i} has an exponential PDF with unit rate. For clarity of notation we drop the subscript *i* in the following proof. It should be taken as implicit.

Proof Assume the last spike was located in bin

*k* and the next spike in bin

*k* +

*L*. If we knew the underlying continuous time conditional intensity function λ(

*t*|

*H*_{t}) and the exact spike time

*t*_{δ} = (

*k* +

*L* − 1)Δ + δ in bin

*k* +

*L* (δ

[0, Δ]) then using the continuous time version of the time rescaling theorem we could write the probability of this event as

which is exponentially distributed in τ. Since we know neither λ(

*t*|

*H*_{t}) nor

*t*_{δ} precisely we must recast τ in terms of what we do know, the discrete bin-wise probabilities

*p*_{k+l}.

The

*p*_{k+l}’s can be written in terms of the underlying continuous process λ(

*t*|

*H*_{t}). Consider any bin

*k* +

*l*. Since discretizaton enforces at most one spike per bin,

*p*_{k+l} does not equal the integral of λ(

*t*|

*H*_{t}) over the bin, but rather the probability (measured from the start of the bin) that the first spike waiting time is less than Δ.

Partitioning the integral in the exponent of

equation 26 into a sum of integrals over each bin allows

*P*(

*t*_{δ}) to be written as.

where we’ve introduced

as a shorthand. By inverting

equation 27 *q*_{k+l} can be written directly in terms of

*p*_{k+l}.

Since we have no information about how λ(

*t*|

*H*_{t}) varies over bin

*k* +

*L* we can pick any functional form, as long as it obeys the constraint that its integral over the bin equals

*q*_{k+L} = − log(1 −

*p*_{k+L}). One choice is

.

^{7} It then follows that

*P*(

*t*_{δ})

*dt* has now been rewritten in terms of what we know, the

*q*_{k+l}’s (implicitly the

*p*_{k+l}’s). We have defined the rescaled time as

(where

) to distinguish it from the rescaled time of the continuous time version of the theorem which directly sums the

*p*_{k+l} and does not require random sampling of the exact spike time δ.

Random sampling of δ

[0, Δ] must respect the choice of

and the fact that we know the spike is in the bin somewhere. For our choice of λ, the probability density of δ conditioned upon there being a spike in the bin is a truncated exponential.

The numerator is simply the event time probability measured from the start of bin

*k* +

*L*. The denominator is a normalization obtained by integrating the numerator between 0 and Δ.

^{8} To draw from this distribution we integrate it to obtain its CDF

set this cdf equal to a uniform random variable

*r* and then solve for δ

This completes the proof.

There are two differences between the discrete and continuous time versions of the theorem. The first, and most fundamental, difference is that

. The latter is only true when Δ is small. Expanding the logarithm of

equation 29 we obtain

To properly rescale the ISIs when Δ is large, all the terms in the Taylor series must be kept. This can be thought of as introducing a correction term (in the brackets) for the finite bin size. Equivalently the approximation 1 −

*p* ≈

*e*^{−p} used in the continuous time version of the proof is not valid for large

*p* (or Δ). The second difference is that we randomly choose an exact spike time

*t*_{δ} = (

*k* +

*L*−1)Δ+δ according to the distribution given in

equation 32. This is done because there is simply no information about where exactly in bin

*k* +

*L* the spike is located and for the rescaled time ξ to be exponentially distributed it must be continuously valued. In the continuous time limit, both of these distinctions vanish.

The hypothesis for testing goodness of fit is now exactly the same as that of the original time rescaling theorem, except that the rescaling is modified to take into account the discretization. Reintroducing the subscript *i* to denote the individual spike times *k*_{i}, the procedure for performing the KS test is simply described.

Procedure for Analytic Correction - Discretize the spike train into bins of width Δ with the spikes in bins {
*k*_{i}} and fit a discrete time statistical model resulting in the spike per bin probabilities *p*_{k}. - Generate a new time series of discrete values
*q*_{k} according to
- For each interspike interval calculate the rescaled ISI ξ
_{i} according to
where δ is a random variable determined by first drawing a uniform random variable *r*_{i} [0, 1] and then calculating
- Make a final transform to the random variables
*y*_{i}
If the discrete time statistical model is accurate the *y*_{i} will be uniformly distributed. Therefore the *y*_{i} can be used to make a KS or Differential KS plot.