Variability and Randomness of the Instantaneous Firing Rate

The apparent stochastic nature of neuronal activity significantly affects the reliability of neuronal coding. To quantify the encountered fluctuations, both in neural data and simulations, the notions of variability and randomness of inter-spike intervals have been proposed and studied. In this article we focus on the concept of the instantaneous firing rate, which is also based on the spike timing. We use several classical statistical models of neuronal activity and we study the corresponding probability distributions of the instantaneous firing rate. To characterize the firing rate variability and randomness under different spiking regimes, we use different indices of statistical dispersion. We find that the relationship between the variability of interspike intervals and the instantaneous firing rate is not straightforward in general. Counter-intuitively, an increase in the randomness (based on entropy) of spike times may either decrease or increase the randomness of instantaneous firing rate, in dependence on the neuronal firing model. Finally, we apply our methods to experimental data, establishing that instantaneous rate analysis can indeed provide additional information about the spiking activity.


INTRODUCTION
One of the primary research areas of computational neuroscience is dedicated to understanding the principles of neuronal coding, i.e., the way information is embedded into neuronal signals. It is generally understood that neurons use brief electrical impulses (called action potentials or spikes) to convey information. The way information is presented in the time sequence of spikes, however, is still a matter of debate (Shadlen and Newsome, 1994;Stein et al., 2005).
A widely accepted answer to the problem is the rate coding hypothesis, which says that the neurons transmit information through the average number of spikes sent along the axon per a certain time window (this is called the mean firing rate). The origin of this theory is credited to Edgar Adrian who found out that the firing rate of the stretch receptor of a frog's muscle changes as a function of stimuli (Adrian, 1926). However, since then, many studies have shown that neurons can encode information without necessarily changing the mean firing rate in response to a stimulus (Perkel and Bullock, 1968;Gerstner and Kistler, 2002;Rigotti et al., 2013;Dettner et al., 2016) prompting the temporal coding hypothesis, which states that the temporal structure of the ISIs is also employed in the embedding of neural information in the spike train (Theunissen and Miller, 1995). Theoretically, information can be encoded in the temporal pattern of the ISIs in an infinite number of ways (Thorpe and Gautrais, 1997;Jacobs et al., 2009;Ainsworth et al., 2012), therefore measures are needed to quantify the features of spiking neuronal signals from different perspectives (Perkel and Bullock, 1968;Victor and Purpura, 1997;Buračas and Albright, 1999;Rieke et al., 1999;Nemenman et al., 2004). A possible way of looking at the role of temporal structures can be through variability coding hypothesis (Perkel and Bullock, 1968) which states that neuronal variability may not be entirely noise, and part of it might reflect the aspects of neural code that is not yet understood; whether it is the variability of the ISIs or of the firing rate. A standard measure for comparing variability in spike trains is through the coefficient of variation (C V ) which is defined as the ratio of standard deviation to the mean of ISIs (Barbieri and Brunel, 2007). Variability of the firing rate is measured by the Fano factor (Ditlevsen and Lansky, 2011;Rajdl and Lansky, 2013;Stevenson, 2016) which is defined as the ratio of variance to the mean number of spikes in a fixed time window.
Another concept that is similar to variability but not equivalent, is the notion of randomness (Kostal et al., 2007b). Variability and randomness both are used as a differentiating measure in the cases where spike trains might seem similar from the rate coding perspective. It is important to distinguish between the two quantities because highly variable data might not be random at all if it consists of relatively predictable values. For example multi-modal data with well separated peaks may have higher variability than uniformly distributed data where the outcomes are the least predictable. Shannon's entropy (Shannon and Weaver, 1949) is widely used to measure randomness (Steuer et al., 2001;McDonnell et al., 2011;Watters and Reeke, 2014), however it is not suitable for continuous distributions. Few other randomness measures based on entropy have been used in neural context recently. In Kostal et al. (2007a) and Kostal et al. (2013) the authors propose a randomness measure for ISIs, creating an alternative to C V ; whereas an entropy based randomness measure of spike counts is introduced in Rajdl et al. (2017), analogous to the Fano factor.
The instantaneous rate is often ambiguously defined as the inverse of a certain inter-spike interval, e.g., of the first complete ISI after stimulus onset, or of a combination of first n intervals, etc. (Bessou et al., 1968;Awiszus, 1988;Lansky et al., 2004). However, statistically consistent definition of the instantaneous rate (Kostal et al., 2018) cannot depend on the "intrinsic" timing given by the particular spike train realization (e.g., the first evoked spike). Instead, it must be evaluated with respect to the "external" time frame, consistently across trials, i.e., asynchronously with respect to individual spike train realizations. In this paper we consider models of spike trains described by renewal processes, and investigate the dependence between the ISI dispersion coefficients and instantaneous rate dispersion coefficients and emphasize that instantaneous rate provides another perspective in the evaluation of neuronal data.
The paper is divided as follows: first, the concept of instantaneous rate is introduced. Next, the concepts of variability and randomness are defined formally. In section 2, we derive the instantaneous rate distribution and the related dispersion measures for a few of the most commonly used models of steadystate neuronal activity. Included models are restricted to the framework of renewal spiking activity for the purpose of this paper. In section 3, we compare the dispersion characteristics of the neuronal models from rate coding and temporal coding perspectives. Upon comparison we find that for some neuronal models, the firing patterns have a different level of randomness in different settings whereas for others, the changed perspective of observation (from rate to temporal) does not affect the randomness of the data. More details on this is available in section 4.

Instantaneous Rate
The class of stochastic processes whose realization consists of a sequence of point events in time is called stochastic point processes (Cox and Miller, 1965). Neuronal spike trains are often described as stochastic point processes, with spikes corresponding to events. There are essentially two ways of describing point processes, either in terms of the number of events occurring in a time window, or in terms of the intervals between these events. Consequently, a spike train can be described either by using a sequence of the occurrence times of its individual spikes, X 1 , X 2 , ... (see Figure 1A) or through the ISIs defined as T i = X i+1 − X i . Generally, the point process is stationary if the underlying joint probability distribution of the numbers of the spikes in k time intervals (Cox and Lewis, 1966, p. 59). In this paper, we consider an important class of stationary point processes, the renewal point processes, in which the length T of consequent ISIs is an independently and identically distributed random variable with the probability distribution function (pdf) f T , T ∼ f T (t). Renewal processes are often used to model the activity of spontaneously active cells (Tuckwell, 1988).
Let the number of spikes that occur in the time interval [0, w] be denoted by N(w). A natural way of calculating the firing rate is to divide the number of spikes N(w) by the time window w. The mean of ISIs, E(T), satisfies the following relationship with the mean of the counting process E[N(w)] (Cox and Lewis, 1966;Rudd and Brown, 1997), where λ is the point process intensity (Cox and Lewis, 1966). For finite w, Equation (1) holds true for renewal processes only if time origin t 0 is arbitrary, i.e., it is not related to the renewal point process realization ( Figure 1B). In this case, t 0 falls into some ISI, say T k . Then the sequence of random variables (W, T k+1 , T k+2 , ...), where W is the time to the first spike from the origin, is not stationary. Moreover, the point process intensity λ is equal to the mean firing rate. The corresponding renewal process is referred to as equilibrium renewal process, as opposed to the ordinary renewal process which starts from an arbitrary spiking event and all the ISIs follow the same renewal pdf (Cox and Miller, 1965).
The instantaneous firing rate is typically defined as the inverse of the ISIs (1/T) (Pauluis and Baker, 2000). However, as proven in FIGURE 1 | Length-biased ISIs and instantaneous rate. (A) An overview of the ISIs when the spikes occur at times X i−1 , X i , X i+1 , .... We assume that the ISIs are independent and identically distributed with pdf f T (t), under steady state conditions. (B) When the observation time t 0 is fixed with respect to some reference time, unrelated to the spike times, the probability of observing a particular ISI is proportional to its length. These "length-biased" intervals (T) are used to define the instantaneous rate R with the property E(R) = 1/E(T). (C) A graphical representation of how the ISI distributions can visually differ from the instantaneous rate distributions, for some well-known ISI models with equal mean firing rate. (D) A summary of the concepts of variability and randomness for the two ways of spike train description that we have considered in this article. Lansky et al. (2004), the mean instantaneous firing rate is higher than or equal to the mean firing rate, with equality if all the ISIs are of the same length. The undesirable inequality in Equation (2) becomes equality once we realize that the "time instant" (at which the instantaneous rate is measured) does not generally coincide with a spike. As shown in Kostal et al. (2018), the instant is chosen with respect to the "external" time frame, across trials, i.e., asynchronously with respect to individual spike train realizations. Consequently, the probabilities of observed ISIs (T) are proportional to their length,T ∼ λtf T (t) ( Figure 1B). The mean inverse of these length-biased ISIs equals to the mean firing rate λ.
The instantaneous rate R = 1/T observed, in this case, is a random variable described by the pdf, For a detailed overview of the derivations, refer to Kostal et al. (2018). An immediate consequence of Equation (3) is, Hence, for the purpose of the derivation of instantaneous rate distributions for different models of steady-state firing ( Figure 1C), we will restrict ourselves to the case of equilibrium renewal process. The variability of instantaneous rate is explored in the next section.

Quantification of Variability and Randomness
The most common method to measure statistical dispersion of ISIs, described by a continuous positive random variable T, is the standard deviation σ (T), defined as The relative dispersion measure quantity is known as coefficient of variation C V (T), where λ = 1/E(T). The coefficient of variation C V (T) is a dimensionless quantity and its value does not depend on the choice of the units of ISIs or on linear scaling; hence it can be used to meaningfully compare the ISI distributions with different means, unlike σ (T) (Softky and Koch, 1993;Doron et al., 2014). From Equation (3), we can write the standard deviation for the instantaneous rate as which leads to the relative dispersion measure, From Equation (5), it follows that σ or C V measure how much the probability distribution is "spread" with respect to the mean value but these quantities do not describe all possible differences between spike trains with equal rate coding characteristics. Spike trains of equal variability may still differ in higher than second statistical moments. Moreover, neither σ nor C V quantifies how random or unpredictable the outcomes are Kostal et al. (2007a).
To quantify randomness as a further describing characteristic of a spike train, entropy based measures like differential entropy (Shannon and Weaver, 1949), h, have been proposed where X is a continuous random variable with pdf f X . However, h(f X ) by itself can not be used as a measure of randomness since it can take both positive and negative values and depends on the scaling of the random variable X. Kostal and Marsalek (2010) proposed the entropy-based dispersion coefficient σ h , Entropy-based dispersion σ h can be interpreted with the help of asymptotic equipartition property (Principe, 2010;Cover and Thomas, 2012), the details of which can be found in Kostal et al. (2013). Analogous to Equation (6), the relative entropy-based measure of dispersion, C h is defined as An immediate consequence of the above equation is that the maximum value of C h is C h = 1, which occurs only in the case of exponential f T .

RESULTS
Among the different point process models of stationary neuronal activity, we have chosen the classically used Gamma, lognormal, inverse Gaussian, shifted exponential, and the mixture of exponential distributions (Tuckwell, 1988). First three distributions are part of the scale family (Casella and Berger, 2002), i.e., their ISI pdfs f T (t; λ), explicitly parameterized by the intensity λ, satisfy the relationship We explore these neuronal models in the following subsections ( Figure 1D). Detailed derivations of the following results are provided in the Supplementary Material for better legibility.

Gamma Distribution
Gamma distribution is frequently used as a descriptor of ISIs in experimental data analysis (Levine, 1991;McKeegan, 2002;Reeke and Coop, 2004;Pouzat and Chaffiol, 2009), its pdf is where Ŵ(z) = ∞ 0 x z−1 e −x dx is the gamma function and a > 0, b > 0 are the parameters. The mean firing rate and the coefficient of variation are, Using Equation (9), the expression for the entropy of the ISI distribution is derived as where Substituting the values from Equations (14) and (15) into Equation (11), gives the dispersion coefficient of randomness, The instantaneous rate distribution f R (r) is obtained through Equation (3) and it follows the inverted gamma distribution, Coefficient of variation C V (R) is evaluated through Equation (8), The differential entropy for f R (r) is, and the dispersion coefficient of randomness is, The dispersion measures C V (T) and C V (R) are related through the following equation, This relationship is illustrated in Figure 2 and we see that Note that the firing intensity λ completely characterizes the exponential distribution and that C V (T) = 1, regardless of the value of λ. The entropy h(f T ) of the exponential distribution is One of the key features of the exponential distribution is that, for a fixed mean firing rate λ, it maximizes the differential entropy among all probability distributions on the real positive half line (Cover and Thomas, 2012). Deriving from Equation (11), The corresponding instantaneous rate distribution, follows from Equation (3), which is the inverse gamma distribution mentioned in Equation (17) with a = 1 and b = λ. For the inverse gamma distribution, the second moment exists only when a > 1. In Figure 3, we can see where gamma distribution becomes exponential and has C h (T) = 1.
FIGURE 2 | Comparison of C V (T) and C V (R) for several standard statistical ISI models. For the gamma distribution, C V (R) → ∞ as C V (T) → 1, since gamma distribution becomes exponential for C V (T) = 1. For lognormal and inverse Gaussian distribution, the relationship between the dispersion coefficients is an identity. For shifted exponential distribution, C V (T) and C V (R) depend on the mean firing rate λ and the refractory period τ . Hence, we vary the value of C V (T) and C V (R) (consequently of λ and τ ) and we see that C V (R) > C V (T) until C V (T) = 0.7715 but then instantaneous rate distribution maintains a higher variability than ISIs. As C V (T) → 1, the shifted exponential becomes exactly exponential and therefore C V (R) → ∞.
FIGURE 3 | Relationship between the dispersion coefficients of randomness for the gamma, lognormal, inverse Gaussian, and shifted exponential distribution of ISIs. Starting from the origin, the arrow indicated to the end corresponds to C V (T) values ranging from 0 to ∞. For the first three distributions, which are part of the scale family and for which we consider λ = 1 for a meaningful comparison, it holds that for C V (T) → 0, we see that C h (T) → 0 for C h (R) → 0. As C V (T) grows, the randomness grows for ISIs and instantaneous rate in the beginning. After a certain point, the randomness starts to decline and as C V (T) → ∞, we have C h (T) → 0 and C h (R) → 0. For the shifted exponential distribution where we have to consider a varying C V (T) and for that a varying mean firing rate λ and refractory period τ , C h (T) increases monotonously as a function of C V (T), whereas C h (R) obtains its maximum for C V (T) = 0.85 and then declines.

Lognormal Distribution
The lognormal distribution represents a common ISI descriptor in experimental data analysis (Levine, 1991;Pouzat and Chaffiol, 2009), even though it is rarely presented as a result of any of the Frontiers in Computational Neuroscience | www.frontiersin.org neuronal models (Bershadskii et al., 2001). The pdf is, where m is the scale parameter and σ > 0 is the shape parameter. The mean firing rate and coefficient of variation are as follows, We compute the differential entropy and the dispersion coefficient of randomness, The instantaneous rate distribution follows the pdf, with The expression for the differential entropy is derived from Equation (11), The dispersion coefficient is evaluated as, For the lognormal distribution, there is a "symmetric" relationship between f T (t) and f R (r) (Kostal et al., 2018), i.e., the shape of the probability distributions of ISI and instantaneous rates are exactly the same for λ = 1. Furthermore, the relationship between C V (T) and C V (R) is that of an identity (Figure 2), and from Equations (29) and (33), we have For lognormal distribution, the randomness and variability of instantaneous rate and ISI are the same, regardless of the perspective, as seen in Figures 2, 3.

Inverse Gaussian Distribution
Inverse Gaussian distribution is often fitted to experimentally observed ISIs (Gerstein and Mandelbrot, 1964;Berger et al., 1990;Levine, 1991;Pouzat and Chaffiol, 2009). The time that a Weiner process with positive drift takes to reach the first passage time is distributed according to the inverse Gaussian distribution. The probability density of inverse Gaussian distribution can be expressed as a function of its mean a > 0 and scale parameter b > 0 From Equation (9), the expression for differential entropy is, where K (1,0) ν (z) is the derivative of the modified Bessel function of the second kind (Abramowitz and Stegun, 1972).
From Equation (8), and from Equation (9), differential entropy is, The expression for the dispersion coefficient of randomness is as follows, Analogous to the lognormal distribution, we observe that the inverse Gaussian distribution also satisfies the "symmetrical" property (Equation 34) and, Results from Figures 2, 3 illustrate this curious property that the randomness and variability of this distribution is equal for ISIs and instantaneous rate perspective.

Distribution Involving a Refractory Period
The refractory period is a state of the neuron, occurring right after a spike, where it is impossible for another spike to be emitted. A shifted exponential distribution function is used as an ISI descriptor for neurons with refractory period τ ( Reeke and Coop, 2004). The probability density function of the shifted exponential function with parameter a > 0 and refractory period τ ≥ 0 is We observe that C V (T) < 1 for a > 0 and τ > 0. The differential entropy for the shifted exponential distribution is evaluated as and substituting these values into Equation (11), we arrive at The pdf of the instantaneous rate is, where Ŵ(s, x) = ∞ x t s−1 e −t dt is the upper incomplete gamma function (Abramowitz and Stegun, 1972). Evaluation of the differential entropy from Equation (9) yields, h(f R ) = − log a 2 1 + aτ − 3(1 + e aτ Ŵ(0, aτ ) + (1 + aτ ) log τ ) 1 + aτ Expression for the dispersion coefficient is derived through Equation (11), For the shifted exponential distribution, we observe that from Equation (48), which leads to, i.e., C V (T) depends on λ and τ . Hence, in order to analyze the relationship between C V (T) and C V (R) we have to vary the values of λ and τ (Figure 2).
Substituting the values from Equations (55) and (56) into Equation (54), we get The dispersion coefficient of randomness C h (R) increases until it attains its maximum value 0.8137 for C V (T) = 0.85, whereas C h (T) keeps increasing monotonously. This behavior can be better understood by looking at the behavior of C h (R) while C V (T) increases in Figure 6C. It is observed that C h (R) attains its maximum value and then declines as C V (T) → 1, whereas C h (T) monotonously increases.

Mixture of Two Exponential Distributions With Refractory Period and Application to Experimental Data
Throughout the years, empirical studies have produced evidence of bimodal or multimodal trends in ISI data (Rodieck et al., 1962;Nakahama et al., 1968;Obeso et al., 2000;Dorval et al., 2008). The underlying assumption is that a neuron might be in one of the several "states, " with each state being characterized by a different ISI distribution (Tuckwell, 1988). A mixture of two distributions is commonly used to modal such data, for example, Bhumbra and Dyball (2004) used a mixture of two lognormal distributions as an ISI descriptor of supraoptic nucleus neurons, and recently gamma-exponential mixture distribution has been used to characterize the ISI distribution in the auditory system (Heil et al., 2007;Neubauer et al., 2009). In this section, we consider the mixture of two exponential distributions (Tuckwell, 1989), which is often used to describe the bursting activity in neurons, where a sequence of short ISIs is dispersed among comparatively longer ISIs. In 1965, Smith and Smith (1965) used the mixture of two exponential distributions to explain the bursting activity in the isolated cortical neurons of an unanesthetized cat. Thomas (1966) used mixed exponential distributions to describe the ISIs in their study of the clustered firing of cortical neurons. Trapani and Nicolson (2011) found that in the lateral line organs of a zebrafish, when the depolarizing currents were blocked, the ISI data of afferent neurons was best described by a mixture of exponential distributions.
The pdf of the mixed exponential distribution with refractory period τ ≥ 0 and mixture components with parameter a > 0, b > 0, a = b is given by FIGURE 4 | Dispersion measures of the mixed exponential distribution with refractory period randomly set to τ = 0.2 and variable weight of components as p ∈ [0, 1] in the direction of the arrows. The rate parameter of one component is fixed at a = 1 whereas the rate parameter b of the second component varies. When p = 0 or p = 1, the distribution is shifted exponential with refractory period τ . The dispersion measure of randomness C h (T) describes C V (T) uniquely when the a = 1, b = 1/2, whereas the reverse is not true (A). This relationship gets complicated with the increasing separation between the rate parameters of the components. As seen in (B-D), the relationship between the other dispersion measures is relatively complicated and non-unique for the given set of parameters.
where p ∈ (0, 1). In this case, The analytical expressions for the C V s and C h s are difficult to obtain, however, they can be calculated numerically for a given set of parameters (Figure 4). The parameter range for this distribution can be vast, we analyze the dispersion coefficients for a few different sets of component rate parameters a and b, given a value of τ and p ∈ [0, 1]. We study the change in randomness and variability as probability variable p increases in the direction of the arrows. The behavior of this model is relatively complicated but for selected cases, such as when a = 1, b = 1/2, or a = 1, b = 1/4, different firing regimes are uniquely described by C h (T) but not by C V (T). Song et al. (2018) model the spike generation in the spontaneously active afferent neurons of the Zebrafish lateral line, as a renewal process. The authors propose that a spike is generated if the neuron is recovered from the refractory period and a synaptic release (excitatory input) from hair cells has arrived and thus ISI T = τ + T E where τ is the absolute refractory period and T E is the time to excitation (we omit the small relative refractory period used in the original paper, FIGURE 5 | Comparison between the dispersion coefficients of variability and randomness of the experimental data for the data sets used in Song et al. (2018). The ISI data has two components: the refractory period and the hair-cell synaptic release time (the time to excitation input). Latter is modeled by a mixture distribution of two exponentials. The color gradient indicates the mean firing rate of the particular data set. The comparison of various dispersion quantities reveal various aspect to differentiate the data sets. As seen in (A), C h (T) classifies the data sets in a similar category from the randomness perspective even if their ISI variability are on a wider scale. The dispersion measure of variability C V (R) also does a better job of differentiating among the data sets in (B) when their randomness dispersion measure would put them in a similar category. In (C), C V (T) and C V (R) reveal separate aspects of the data sets. Overall C V (R) helps differentiate among the data sets with similar C V (T) values. Finally in (D) we can see that C h (R) further differentiates the data sets with equal C h (T). The data sets which might have similar randomness on the ISI scale, can be differentiated on the basis of their instantaneous rate randomness. All of these figures support our assertion that dispersion measures of instantaneous rate provide additional information that can be helpful in distinguishing the data sets.
for computational simplicity). It is demonstrated in the paper that the mixture of exponential distributions used to model the hair cell synaptic release time T E , yields the best fit for the ISI data. The pdf of the ISI, in this case is given by Equation (58). We calculated C h and C V for ISIs and the instantaneous rate of each data set fitted with a combination of absolute refractory period and mixture of exponential distributions. For the given data sets, C V (T) and C V (R) offer more information than their randomness counterparts C h (T) ( Figure 5A) and C h (R) (Figure 5B) respectively. When it comes to the dispersion measures of variability, as seen in Figure 5C, in some situations C V (R) offers additional information to distinguish further nuances in the data. The data sets which might seem of similar ISI variability, differ when it comes to the variability of their instantaneous rate. Similar results follow for the randomness (Figure 5D). This example supports our assertion that the dispersion measures based on the instantaneous rate can provide additional information to help differentiate among the data sets.

DISCUSSION
We studied the spiking activity described by the renewal processes from two perspectives, the temporal point of view (in terms of ISIs) and the frequency point of view (in terms of instantaneous rate). We found that for a given spike train the temporal characteristics and the instantaneous rate characteristics, can either follow the same trend or opposite trends. This is due to the fact that the instantaneous rate distribution is obtained from length-biased sampling of ISIs (Equation 3). Spike trains described by non-renewal processes have been studied widely (Eden and Kass, 2016) but are beyond the scope of this paper. For the special case of serially correlated ISIs, the results of our analysis apply to marginal distributions and therefore remain unchanged (Kostal and Lansky, 2006). 1. We observe several different relationships between the variability from temporal and instantaneous rate perspectives (Figure 2). For gamma distribution, the variability of instantaneous rate is higher than the variability in ISIs whereas for lognormal and inverse Gaussian distribution it stays the same, C V (T) = C V (R). On the other hand, for shifted exponential distribution C V (R) < C V (T) until C V (T) = 0.7715 but after that C V (R) increases rapidly compared to C V (T). 2. In the case of gamma distribution, both the randomness measures C h (T) and C h (R) decline eventually as a function of C V (T) but they become constant as a function of C V (R) (Figures 6A-D). The randomness for gamma distribution differs in each particular case, and for small values of C V (T), C h (R) < C h (T). Lognormal and inverse Gaussian distributions attain their maximum for C V (T) values close to 1 and then their randomness keeps declining (Figures 6A,C). For these two distributions, C h (T) = C h (R) whether it is mapped against C V (T) or C V (R). For shifted exponential distribution, C h (T) increases as a function of C V (T) ( Figure 6A) whereas C h (R) attains its maximum for C V (T) = 0.85 and then declines fast as C V → 1 (Figure 6C). 3. We studied the case of mixed exponential distribution with a refractory period. Although the theoretical analysis is complicated, we use the experimental data obtained from Song et al. (2018) to inspect the temporal and instantaneous rate perspectives. The dispersion measures for instantaneous rate provide a novel outlook on the data, different from the one provided by the dispersion measures for the ISIs (Figure 6). In cases like these, the instantaneous rate may be helpful in distinguishing further nuances in the data.

AUTHOR CONTRIBUTIONS
LK proposed the concept of the article. RT and LK wrote the manuscript. RT performed the calculations, numerical simulations, and analysis. All authors contributed to the article and approved the submitted version.

FUNDING
This research was supported by the Czech Science Foundation project 20-10251S.