Fano Factor: A Potentially Useful Information

1. Introduction

The frequency by which neurons generate spikes (action potentials) is commonly considered as a basic form of information transfer within the neuronal system (Adrian and Zotterman, 1926; Perkel and Bullock, 1968). This “frequency code” is often quantified by the number of spikes within an appropriately set time window. The length of the window is limited by the requirement of stable conditions on one side, and by aiming at reproducible results on the other side. The natural conditions typically vary rapidly, but even if kept constant by an experimenter, there are other reasons, like spiking adaptation (Benda and Herz, 2003), which restrict the duration of the observation time window. All these constrains create difficulties for the statistical inference based on the firing rate, and many sophisticated methods to overcome them have been developed (Nawrot et al., 1999; Dayan and Abbott, 2001; Cunningham et al., 2009; Benedetto et al., 2015; Kostal et al., 2018; Tomar, 2019).

As mentioned, the step from the observed times of spikes to the firing rate quantification requires a proper selection of the statistical procedure. The data available from experiments are usually formed by repeated spike trains (repeated trials) recorded from a single neuron exposed to the identical conditions (reaction to a stimulus). It aims to reflect activity of a neural network of a statistically homogeneous population of neurons. In other words, multiple trials are only an experimental counterpart of real-time multi-neuronal activity. Experiments of multi-unit recordings are also available, however, there is currently no technical tool to simultaneously record all spike trains coming to a target neuron. Independently, whether data come from the repeated experimental trials or from a simultaneous multi-unit neuronal recording, in both cases additional assumptions are often accepted, namely, that the data are independent and statistically identical replicas of the same random variable (which is also assumed in this paper). Then, the common definition of the firing rate is the mean number of spikes per unit time window.

It has been noticed for a long time, and under various experimental conditions, that the firing rate (however defined) does not completely characterize neuronal activity, as demonstrated by many studies going beyond the rate coding concept (Shadlen and Newsome, 1994; Rieke et al., 1999; Olypher et al., 2003; Stein et al., 2005; Kostal et al., 2007; Christodoulou and Cleanthous, 2011; Rajdl et al., 2017). The “neural codes” based on rate-independent components are generally called temporal codes. Probably the simplest one going beyond the mean number of spikes is the variability of that number.

Similar to the firing rate, there are several possibilities on how to quantify the variability of spike trains or in generally non-biological systems (Lindner et al., 2004, 2005; Bravi et al., 2011; Kostal et al., 2013; Aoki et al., 2016; Bowden, 2017; Peterson and Heil, 2018; Ilan, 2020). The variance of the spike count in the observation window would be the typical first candidate. However, large numbers usually have large variance (the second moment is quadratic in the scaling parameter), hence some “relative” dispersion measures are often more desirable. Based on properties of Poisson distribution, as the Poisson process is some kind of template for any series of uniform events appearing in time, the index of dispersion was introduced (Cox and Lewis, 1966). It relates the variance of the number of spikes in a time window to its mean. In Fano (1947) a metric for directly unmeasurable quanta was proposed. It was based on the Poissonian character of the quanta appearance and was also represented as the ratio between variance and mean (Fullagar et al., 2017). It became known as the Fano factor. Despite originally being unrelated to the time evolution and not at all admitting data with a non-Poissonian character, due to the formal equivalence with index of dispersion, the term Fano factor prevailed regardless of the difference in aims. The term Fano factor will therefore be used throughout the paper.

For the Poisson distribution the variance-to-mean ratio is equal to one. This fact has fascinated neuroscientists for many decades and the Fano factor has been evaluated and presented over a wide range of experimental conditions, different types of neurons, and species of experimental animals. The number of references to papers in which it was evaluated would be practically endless. Most of the theoretical models of neurons are oriented toward the calculation of interspike interval (ISI) distribution and directly assume a renewal character of the firing. Many of these studies present coefficient of variation of the interspike intervals. Considering that the squared coefficient of variation is equal to the Fano factor (over an infinite window), experimental studies and also countless numbers of the theoretical papers present the Fano factor as a property of the models investigated in them.

As mentioned above, the mean number of spikes in the observation window is the most common metric for the firing rate. Thus, relating the variance to this quantity may induce the wrong feeling that the measure (index of dispersion, Fano factor) is firing rate independent (Figure 1). Such independence is valid only asymptotically for exceptionally fast firing and exceptionally slow firing, in both cases with respect to the time window over which the counts of spikes are investigated. Neglecting rate changes while studying the Fano factor can lead to incorrect, distorted results, mainly when comparing the Fano factor based on various sets of spike trains from various experiments. A change of the spiking rate changes the estimated values of the Fano factor even if the true values remain constant. This property of the Fano factor is known and there are some ways to avoid this problem, e.g., using the operational time (Nawrot et al., 2008) or the mean matching method (Churchland et al., 2010). However, it seems to us that this issue is still often neglected and has not yet been sufficiently described and explained. For examples of an insufficient Fano factor analysis see the Discussion in Nawrot et al. (2008). Our aim is to describe the dependence of the Fano factor on the intensity in more detail, to show possible problems, and to present and test a method on how to avoid them. The main general solution is to use the Fano factor in the operational time, for the firing rate equal to one. The above-mentioned mean matching method can be suitable too, however, it is usable only in specific situations and is thus not investigated here.

The paper is organized as follows. The relevant theory is summarized in the first part of the paper. We introduce renewal processes and Markov renewal processes, which we use to model the spike trains. Then, the influence of the firing rate on the Fano factor is described and illustrated and a rate-independent variant of the Fano factor is defined. Finally, a method for comparison of variability of two or more data sets with different firing rates is proposed and evaluated.

2. Model of Spike Trains

The Fano factor deals with the number of events (spikes) within an observation time window and thus we have to specify the basic properties of such a counting process. The equilibrium counting process N(w) describes the number of spikes in an interval (0, w], w > 0, where time zero is “randomly positioned” with respect to the sequence of spikes. The dynamics of the process is captured by its intensity function, λ(w), defined by equation

which reflects the rate for the occurrence of a new event for a short time window [w, w + dw). Under the assumption of stationarity, which means that all the statistical properties of N(w) are time invariant, it holds

where E(T) is the mean of ISIs, which due to the stationarity of the sequence of ISIs is constant, E(T) = μ. Combining (1) and (2) we can see,

Unfortunately, formulas analogous to (2) relating higher moments of the counting process to the higher moments of ISIs do not exist. The counting process N(w) is fully determined by the probabilities p_n(w) = P(N(w) = n)). If N(w) has Poisson distribution with parameter λw, the process is called the Poisson process. For other commonly investigated point processes there is no analytical expression for p_n(w).

There is the basic duality relationship between the counting process, N(w), and intervals between events (ISIs),

where T₀ is time of the first spike and T_n−1 is the (n − 1)th complete ISI. The spiking process is called the equilibrium renewal process under the additional assumption that the ISIs are independent and identically distributed continuous and positive random variables, denoted by T, with the probability density function (pdf) f(t). The variants usually employed for pdf f(t) will be summarized later. The term “equilibrium” again specifies that the time zero is unrelated to the sequence of spikes.

Besides the intensity λ, a basic and very often used characteristics of T (ISIs) is the coefficient of variation, defined as

where Var(T) is the variance of T. The main advantage of CV of ISIs, as a measure of spike train variability, is that CV is dimensionless and that ISI probability distributions with different means can often be compared meaningfully. More precisely, CV does not depend on λ if the ISI pdf belongs to the “scale (or rate) family,” i.e., if the pdf f(t; λ) (explicitly parameterized by the rate λ) satisfies f(t; λ) = λf(λt; 1) for any λ > 0 (Casella and Berger, 2002). For example, the exponential pdf with fixed refractory period, discussed later in the paper, does not belong to the scale family. Similarly, the value of CV is not changed if the time axis is scaled linearly (e.g., if the refractory period is also rescaled as in our example in section 4).

Renewal processes are standardly used to model neuronal spike trains (Shinomoto et al., 2005; Nawrot et al., 2008; Shimokawa et al., 2010; Omi and Shinomoto, 2011; Ostojic, 2011; Fisch et al., 2012; Pipa et al., 2013; Koyama and Kostal, 2014; Lansky et al., 2016), however, sometimes they might be seen as too simple. We thus also use a generalization of renewal processes–Markov renewal processes (MRPs). In this case the ISIs are not identically distributed, but every ISI has one of the n ∈ ℕ probability distributions with pdfs f_i(x), i = 1, …, n. The pdfs of the individual ISIs are determined by the states of a Markov process, with a transition matrix P = (p_i,j), 0 ≤ p_i,j ≤ 1, i, j = 1, …, n, expressing the probabilities that after an ISI with pdf f_i(x) there will be an ISI with pdf f_j(x). Such a process is a combination of a Markov chain with n states and a renewal process. Analogously to renewal processes, we assume MRPs in equilibrium. For more detailed description of MRPs see (Çinlar, 1969; Cox and Isham, 1980).

For simplicity, we focus only on a variant with n = 2 and p_1,2 = p_2,1 (the probability of transition from the state one to the state two is the same as the probability of reverse transition), thus we have two random variables with pdfs f₁(x) and f₂(x) (with means μ₁ and μ₂ and coefficients of variation CV₁ and CV₂) and a transition matrix

where 0 ≤ p ≤ 1 is the probability that the state will be changed. The intensity of such a process is

which yields $\lambda /2<{\mu }_1^{-1}\le \lambda \le {\mu }_2^{-1}$ (for μ₂ ≤ μ₁). If the difference of μ₁ and μ₂ is large and p is small, we obtain a bursting process. On the contrary, if p = 1, the process always changes the states (distributions of ISIs) and is called an alternating renewal process (ARP). A basic example of MRPs (resp. ARPs) is the Markov (resp. alternating) Poisson process (MPP, APP), when the ISIs have exponential distributions $f_1(x)={\mu }_1^{-1}exp\left\{-x/{\mu }_1\right\}$ and $f_2(x)={\mu }_2^{-1}exp\left\{-x/{\mu }_2\right\},{\mu }_1,{\mu }_2>0$.

Finally, let us mention a special and interesting case of equilibrium renewal processes situation for CV → 0. We obtain a sequence of equidistant points (times of spikes) with a random origin. All the ISIs are thus constant with length 1/λ and the only source of variability is the time to the first spike, T₀, which is distributed uniformly. This model is called the (equilibrium) pacemaker (PM) and it represents the limit case of renewal processes from the point of view of (low) variability.

3. Fano Factor

Fano factor is a measure of variability of a counting process defined as

thus, as the variance to mean ratio of the number of spikes in a time window of a length w. As mentioned, this quantity was originally called the index of dispersion for count (Cox and Lewis, 1966). Fano factor defined by (8) is a function of w, however, the same term is often used to directly denote the limit

which removes the dependence on w.

For renewal processes, basic properties of Fano factor are (Cox, 1962)

Equation (11) has such an effect that there are cases when CV² is also called Fano factor (Shuai et al., 2002). The function F(w) cannot be expressed analytically for most of the probability distributions of T, however there are some formulas that can be solved numerically, e.g., Rajdl and Lansky (2014),

where ${L}$ and ${{L}}^{-1}$ denote the Laplace transform and its inversion. For CV → 0 (the equilibrium pacemaker process), we obtain

where τ = ⌊λw⌋. The limit of formula (13) for w → ∞ is (as expected) zero.

For the Markov renewal processes, analytical calculation of F(w) is even more complicated. At least the limit of F(w) for w → ∞ can be derived in a closed form, using the results from Ball and Milne (2005),

The limit for w → 0⁺ equals one, as for the renewal processes. For MPPs relationship (14) reduces to

and for APPs to

The parameter p thus highly influences F, and by decreasing p we can increase F arbitrarily. Also, an alternating Poisson process is always more variable than a simple Poisson process; from Equation (16) we can see 1 ≤ F < 2, with F = 1 only for μ₁ = μ₂ (standard Poisson process). Equations (7) and (15) yield that to obtain a MPP with given λ and F (for a fixed p), the individual means need to be

Estimates of the Fano factor, F(w), and the intensity λ are standardly calculated based on observed numbers of spikes N₁(w), N₂(w), …, N_n(w) in a time window (0, w], w > 0, where n is the number of repeated trials or number of simultaneously recorded neurons (see Figure 2A). The standard estimators are

Formula (20) is clearly an estimator of F(w), however, it is used to estimate also directly the limit F.

The mean square error (MSE), the common quantity to measure the quality of the estimator, can be decomposed to bias and variance,

The bias and variance of the Fano factor estimator can be approximated as

Formula (22) shows that the main part of bias of the estimator is unsurprisingly created by the difference F(w) − F.

Besides MSE, a standard measure of error of an estimator is the mean absolute error,

which is a more natural and better interpretable error measure than MSE, however, it is more difficult to analyze it theoretically.

4. Influence of Intensity on Fano Factor

The purpose of the Fano factor evaluation is to measure the variability without the effect of the intensity. However, only the limit value F is generally intensity independent. To show this fact, let us consider two identical equilibrium renewal processes N₁(w), N₂(w), which differ only in the values of intensities λ₁ > λ₂ > 0. For this process, it clearly holds

which directly yields

The values of the Fano factor for the processes N₁(w) and N₂(w) are thus generally not equal for all w > 0, a change of the intensity changes the value of F(w) even if all other characteristics are constant.

For an illustration of these facts we show the dependence of F(w) on the intensity λ for renewal processes with several very often used probability distributions of ISIs, gamma, inverse Gaussian, and shifted exponential. These distributions can be parameterized using λ and F, their probability density functions are then

for the gamma distribution, where Γ(x), x > 0, is the gamma function, and

for the inverse Gaussian distribution. The most basic distribution of ISIs is the exponential distribution,

yielding the sequences of spikes corresponding to a Poisson process. However, as for this distribution F(w) = 1 for any w > 0, there is no dependence of Fano factor on the intensity. Situation changes when we extend the model for a refractory period (RP) of a length r ≥ 0. Refractory period is a time period after each spike when the probability of occurrence of another spike is zero. ISIs then correspond to a random variable T = r + T′, where T′ has exponential distribution. If we want such a process to have intensity λ and Fano factor F, F ≤ 1, its density has to be

and zero otherwise.

To calculate the values of F(w) we numerically solved formula (12) and confirmed the results by Monte-Carlo method using (20). We also show the values of F(w) for MPPs with p = 1 (APP) and p = 0.1 (a bursting process). In these cases, we calculate the means μ₁ and μ₂ using Equations (17) and (18) and the values of F(w) using Monte Carlo simulations. Note that the results for these processes are shown only for F = 1.5, as F cannot be lower than one. The limit case with zero Fano factor (the pacemaker model) is shown too.

The results are shown in Figure 3, and we can see that the dependence on the intensity is high, yielding various values of F(w) without changing w or F. As already mentioned, the dependence on λ is the same as the dependence on w, thus we obtained known curves showing the dependence of the Fano factor on the length of the observation window. Note that this holds assuming that the refractory period, included in T, scales with the firing rate. If the refractory period is fixed, not influenced by the changes of intensity, the interchangeability of w and λ does not hold. There are several interesting features, we would like to point out. Firstly, it is the non-monotonicity of the curve for inverse Gaussian distribution with F = 1.5. As shown and explained in Rajdl and Lansky (2014), this holds for any distribution with f(0) = 0 and F > 1. Secondly, we can see that F(w) for gamma distribution with F = 1.5 converges to unit value for λ → 0 very slowly. It is caused by the fact that probability of occurrence of two spikes in a very short window is in this case not negligible, which results from the specific shape of the gamma distribution with F > 1 (the probability of occurrence of a spike near zero is relative large). Finally, for the peacemaker model, the dependence is non-monotonic, with a cyclic-like behavior. Fano factor is zero in integer multiples of the length of the window w, when we know exactly how many spikes will be observed.

5. Comparison of Fano Factor in Experiments With Various Spike Intensity

It follows from the previous section that the values of F(w) cannot be studied separately without considering the values of intensity λ. The natural way how to deal with this problem is to compare values of Fano factor only under the condition that the intensity is fixed. Generally, it is suitable to transform time so that the intensity equals one, before calculation (estimation) of F(w). Such a time (window) is called operational time (Nawrot et al., 2008) and can be expressed as w_o = wλ. Two values of Fano factor can be then compared if calculated for the same w_o (see Figure 2B). To avoid prior explicit time transformation, we can define operational Fano factor as

where w_o is the operational time (time expressed in multiplies of mean ISIs), see Equation (1), but F(w_o/λ) is calculated on the original time scale. Characteristics F_o(w_o) is then a variability measure fully independent from the intensity (for equilibrium renewal processes) and it is thus directly comparable even for different experiments.

For estimation of quantity (31) we need to also estimate the intensity. It can be easily done by estimator (19), which yields the estimator of operational Fano factor in form

Measuring Fano factor in operational time is especially important if we want to compare its values estimated based on data from various experiments with various conditions. Nevertheless, there are also some issues complicating this method. Firstly, while estimating the operational Fano factor, we work only with estimated intensities. Estimation of the intensity is a simple task, but it still brings some additional uncertainty (variance) to the value of estimator (32), which could reduce its accuracy. Secondly, we need to determine and use only the greatest common operational time ${\overline{w}}_o$ available for all the data (experiments). Suppose that we have m > 0 experiments with observation windows w_i and estimated intensities ${\widehat{\lambda }}_i$, i = 1, …, m, then

see Figure 2B. This implies that some data (spiking times) might not be used, which could seem contra productive. To clarify the severity of these problems and positives of the operational Fano factor while comparing values from different experiments, we perform various evaluations based on Monte Carlo simulations.

In the following we suppose that we have measured the responses, times of spikes fired in a window (0, w], of a neuron under two different stimuli, and we are interested in the change of variability in the responses, measured by Fano factor. Thus, we have two sets of n > 0 spike trains, first set with intensity λ₁ and Fano factor F₁ and second with intensity λ₂ and Fano factor F₂. We aim to find out, how reliable is the ratio r(F) = F₂/F₁ estimated using the standard estimator (20) and using the estimator of the operational Fano factor (32). Let us denote the estimated ratio for the standard Fano factor as $r(\hat{F})$ and for the operational Fano factor as $r({\hat{F}}_o)$.

We study this situation based on data generated using several models of spike trains with various parameters. Specifically, we use equilibrium renewal processes with gamma and inverse Gaussian distribution, MPP with p = 0.1 and MPP with p = 1 (APP). Without the loss of generality, we fix λ₁ equal to one and vary λ₂ in the interval (0, 5]. For the length of the window in operational time we use values 1, 5 and 10 and always generate n = 50 spike trains. For F₁ and F₂ we use values 0.5 and 1.5 in all combinations, thus we assume situations when Fano factor increases, decreases and remains constant. Every situation was repeated (generated) 2,000 times.

While using the operational Fano factor, the observation window in the set of spike trains with the larger value of $w\widehat{\lambda }$ is shortened. It is thus possible not the use the window which starts at zero, but to shift it within the original window. For example, besides starting at zero, a reasonable possibility could be shifting it to the middle or select its position randomly. This decision depends mainly on where the most likely information is that we want to capture. Nevertheless, while using the equilibrium processes as the spike trains model, the location of the window does not matter. Therefore, in the following examples we use zero as the beginning of the sampling window. Another possibility, how to deal with this question, is to estimate Fano factor repeatedly while shifting the shortened window gradually in the original one and use average as the final estimate. This method also removes the disadvantage of the operational Fano factor, regarding omitting some time of spikes after shortening the window. We evaluate performance of this method too, specifically, we shift the window gradually so that all the spikes are used, but as few as possible spikes are used twice. We denote the corresponding estimator ${\hat{F}}_{õ}$.

Firstly, in Figure 4, there are shown the median values (across the 2,000 repetitions) of $r(\hat{F})$ and $r({\hat{F}}_o)$ in situation when F₁ = F₂. In this scenario, the true ratio of the Fano factor values equals one. While changing intensity λ₂, the standard Fano factor shows systematically non-unit ratios, thus a false change of variability. The problems are larger mostly for smaller w, when the dependence of F(w) on λ is stronger. Further, the ratios do not tend to unit value for increasing λ₂, as F₁(w) is generally different from F₁ [to which converges F₂(w)]. Note that the initial slight non-monotonicity of the curve for the inverse Gaussian distribution in situation with F₁ = F₂ = 1.5 and w = 1 is caused by the non-monotonicity of the corresponding curve in Figure 3. Contrary to these inconvenient properties of standard Fano factor, the ratios of the operational Fano factor are always (near) one, the dependence on the intensity is removed. It holds even for MPP with p = 0.1, where the dependence of standard Fano factor on the intensity is very high.

When the real values of Fano factor remain constant (F₁ = F₂), operational Fano factor fully removes the potentially false change of estimated r(F), as expected. However, it does not immediately imply, that the operational Fano factor gives better information about r(F), e.g., the estimator can have larger variance (as explained earlier) and thus possibly also larger MSE or MAE. To explore this possibility we show, in Figure 5, by how much is MAE of the standard estimator larger than MAE of the estimator of the operational Fano factor [assuming one as the real value of r(F)]. We also show the difference from the improved operational Fano factor $r({\hat{F}}_{õ})$. We can see that the operational Fano factor almost always decreases MAE, the few opposite cases are likely just random events. Moreover, the estimator $r({\hat{F}}_{õ})$ further decreases MAE significantly.

In Figure 6, there are shown some situations where the values F₁ and F₂ in the individual experiments differ. In this case, it is more complicated to determine the optimal results. Of course, ideally, we would like to obtain value F₂/F₁ independently of λ₂ as the ratio of the Fano factor estimators, but it is clearly not possible (it would require w going to infinity). The best result we can hope for is to obtain constant ratios of Fano factor estimates (independent of λ₂), positively correlated to F₂/F₁. As we can see, the standard Fano factor does not fulfill this requirement, the differences highly depend on λ₂ (mainly for small w). The behavior of the operational Fano Factor is better, as the ratios are intensity independent at least for λ₂ > 1. However, for λ₂ < 1 they change highly even for this variant of Fano factor. It is caused by the fact that decreasing of intensity decreases usable length of the window in operational time, which causes convergence of Fano factor to unit value. The difference from the Figure 4, where shortening of the operational window causes no problems, is that in Figure 4 the bias of the Fano factor estimates (for λ₁ and λ₂) due to the short operational window is the same for the both estimates and therefore cancels out in their ratio. The transformation to the operational time removes any bias of the Fano factor estimate caused by differences in intensities, but there is another type of bias due to the limited usage of spikes, which is enhanced when the window is short.

6. Real Data Example

Similar to the previous section, we study the change of variability between two sets of spike trains but now recorded from real neurons under different experimental conditions. Our aim cannot be to show that information about variability deduced from the operational Fano factor is more accurate than that from the standard method because, in contrast to the simulated data, we do not know the truth. We also do not intend to deduce any specific conclusions about the data, our aim now is only to show that these two methods can yield very different results. The data sets were obtained from www.crcns.org (Kohn and Smith, 2016), where can be found details about the experiment. Briefly, the data consist of spike trains measured in visual cortex of three adult Macaca fascicularis monkeys under three different visual stimuli. While measuring the neuronal activity, the monkeys (anesthetized and paralyzed) were watching three different movies (white noise, a natural scene, and sinusoidal gratings), which yielded spike trains of length 30 s from dozens of neurons. Each experimental setting and measurement were repeated identically 120 times. For our purpose, a much smaller data set is sufficient, we use only a specific part of the data. We compare Fano factor in two of the three experiments, monkeys watching the white noise (experiment 1) and the natural scene (experiment 2). We also use only data of one monkey from time interval 5–10 s as it is sufficient duration for our purpose. This part of data was selected rather randomly as exploring other parts does not yield fundamentally different results. In this data, there are 74 neurons each with 120 repeated spike trains. Based on them we calculate the change of Fano factor in various windows of length w and the change using the operational Fano factor concept. The operational Fano factor is calculated in window ${\overline{w}}_o$, using formula (33). Thus, we first estimate the intensities in window w, ${\widehat{\lambda }}_1$, and ${\widehat{\lambda }}_2$, and then shorten the window in the experiment with larger intensity according to the ratio of the intensities. These calculations are performed for each neuron separately—the spike trains from experiments 1 are compared to the spike trains from experiment 2 for the same neuron. An example of the data is shown in Figures 7A,B. Times w and ${\overline{w}}_o$ illustrate shortening of observation window (w) in the spike trains with larger intensity, to obtain the greatest common operational time (window) ${\overline{w}}_o$.

This way we obtain (for every w and neuron) two estimates of standard Fano factor, ${\hat{F}}_1$ and ${\hat{F}}_2$, and operational Fano factor, ${\hat{F}}_{o1}$ and ${\hat{F}}_{o2}$. As in the previous section, we measure the change of the estimates as ratios, $r(\hat{F})={\hat{F}}_2/{\hat{F}}_1$ and $r({\hat{F}}_o)={\hat{F}}_{o2}/{\hat{F}}_{o1}$. Note that median values of Fano factors and intensities estimated in the whole window (w = 5s) are $m({\hat{F}}_1)=6.62,m({\hat{F}}_2)=4.05,m({\widehat{\lambda }}_1)=3.48\text{spikes}/s$, and $m({\widehat{\lambda }}_2)=4.77\text{spikes}/s$. The values of F are relatively high, which is given partially by bursting character and partially by some inhomogeneous behavior of the spike trains in time, resulting from the nature of the experiments (mainly for the movie scenario). This violates the theoretical assumptions, however, as we study the data mostly in short time intervals, the influence should not be serious.

The relation of $r({\hat{F}}_o)$ to $r(\hat{F})$ is shown in Figure 7C. The value d is calculated as the mean absolute difference between $r(\hat{F})$ and $r({\hat{F}}_o)$. We can see that the differences of the results obtained by these two methods can be relatively large, there are even frequent cases when one method shows increase of variability and the other one a decrease [the largest difference is for $r(\hat{F})=1.88$ and $r({\hat{F}}_o)=0.78$]. Of course, for small w the differences are small (w = 0.002 and w = 0.005), but they increase with increasing w and later decrease again. For w = 5 the difference d is practically the same as for w = 0.005 (although the distributions of the underlying data are different). It corresponds to the theoretical results, for both small and large w the influence of intensity on Fano factor decreases and thus the application of the operational Fano factor in these situations does not change the results much. We can also see that for ${\widehat{\lambda }}_1<{\widehat{\lambda }}_2$ it is more likely that $r(\widehat{F_o})<r(\hat{F})$ and vice versa. This is not a general rule, in this case, it is caused by the fact that the spike trains are mostly highly variable (F > 1), which leads to this behavior. The reason is that we are shortening the window in the second experiment and then, from Equation (10), the value of Fano factor mostly decreases nearer to one, decreasing also the ratio $r(\widehat{F_o})$.

7. Discussion and Conclusions

The dependency of the Fano factor on firing intensity has been investigated in the previous sections. We studied the effect of the firing rate on the sampling Fano factor F(w). This approach should not be confused with the study of the objective dependency of the Fano factor, F, on the firing intensity, which is described in some experimental studies and in many theoretical models. If an increase in the firing rate is not accompanied by the same increase in the variance of the number of spikes, then the Fano factor F changes with the firing rate. For theoretical treatment of this problem, see Koyama (2015), Stevenson (2016), and Charles et al. (2018).

The Fano factor is commonly seen as an analogy to CV based on the numbers of spikes instead of ISIs. They both create two basic measures of variability in series of events. However, the analogy is not accurate. The reason is that the idea of removing the influence of the rate does not work for F(w) as for CV as the random variable N(w) depends on linear change of rate (time). The definition (8) is based on a different justification, namely, it is the variance of a counting process related to the variance of Poisson process with the same intensity. For Poisson process the mean count equals its variance and it implies Equation (8). Because of the importance of Poisson process, F(w) is just a relative variability of the investigated process with respect to the Poisson process. The reason F(w) became popular as variability measure is its relation to CV for renewal processes (11). However, there is an important difference from CV, the additional parameter—the length of the observation window w, which (often highly, mainly for lower w) influences the value of F(w), whereas (theoretical) CV remains constant. Of course, in practice, estimation of CV is also performed for a limited number of ISIs (in other words, in a fixed observation window), leading to similar problems, but of a different nature (caused by censoring of ISIs) (Pawlas and Lansky, 2011; Rajdl and Lansky, 2014). The aim of this paper is to stress that the Fano factor does not parallel the simple notion of CV of ISIs for the case of counts, as often implicitly assumed. We show that even for simple and commonly used renewal models the dependence on λ must be accounted for (i.e., by the method of operational time).

If the Fano factor is employed to measure the variability of spike trains (or other events occurring in time), the effect of the intensity has to be removed. This is not an unknown fact, but it is almost always overlooked or (probably) not entirely understood. We have explored and explained this issue in detail and showed that Fano factor needs to be used carefully. First, we have theoretically explained the influence of intensity on the values of the Fano factor, and illustrated it using several models of neural spike trains. Next, we described how to avoid potential problems - by estimating Fano factor in operational time and defined the intensity-independent approach (the operational Fano factor). Finally, we showed the benefits of the operational Fano factor while comparing the variability in two experiments with different intensities. For the theoretical results and studies based on Monte-Carlo simulations, we used renewal processes and their generalization (Markov renewal processes) as the spike train models.

It is important to emphasize that the theoretical part of the paper was based on specific models of the spike trains, mostly renewal processes. In reality the behavior of neuronal spike trains can be of course much more complex. Let us remind here the main simplifications made when modeling spike trains using the renewal processes, as their violation can cause misinterpretation of the achieved results, for a more detailed discussion see Nawrot et al. (2008). Firstly, the stationarity in time is assumed. This condition can be sufficiently fulfilled in spontaneous activity, but seriously violated in a response to a stimulus. In such a situation, estimation of time dependent intensity and following transformation to operational time can be performed as in Nawrot et al. (2008). However, when the time dependence of intensity is not too significant and the observation window is relatively small (as it mostly is while estimating Fano factor), assuming stationarity in time can be sufficiently accurate. The second possible violation of the theoretical assumptions is non-stability across trials, e.g., when the firing rate in consecutive trials is not constant. This is a troublesome problem, as it is difficult to distinguish between the intrinsic changes and systematic ones. An indicator of this issue can be joint analysis of both F and CV, as in such a situation F tends to be larger than CV² (Nawrot, 2010). Finally, the renewal assumption is often found violated in real data and it has a large effect on the dependency of the F estimate on the length of the observation window, as shown in Chacron et al. (2001), Farkhooi et al. (2009), and Avila-Akerberg and Chacron (2011).

Additionally, to these three potentially critical points we assumed that scaling of time is equivalent to the scaling of intensity, in other words that the ISI distribution belongs to the scale family. It is also a simplification, however, the majority of frequently applied ISI distributions (Gamma and Inverse Gaussian) belong to this family. A warning example when this condition is violated can be consideration of the absolute refractory period, expressed in real (fixed) time. It would imply that the intensity and the observation window are not fully interchangeable and estimating both standard and operational Fano factors, these issues should be kept in mind.

Under the mentioned assumptions, the results can be summarized in the following points.

• Values of Fano factor calculated in a finite window depend on the intensity, even if other characteristics of the spike trains are fixed. Intensity independent (for renewal processes) is only the limit value F, the values F(w) are influenced by λ, mainly for small w.

• the influence of the intensity is caused by the fact, that changes of intensity are equivalent to changing the length of the observation window with fixed intensity. Note that this holds under the definition of intensity λ = 1/E(T) and T from the scale family of distributions. However, for example, in a model with absolute refractory period, the interchangeability of w and λ does not hold.

• to remove the influence of the intensity, it is suitable to calculate Fano factor in time linearly transformed so that λ = 1 (operational time) and we show how to calculate an operational Fano factor.

• the influence of intensity has to be taken into account mainly while comparing results from two (or more) experiments. If the intensity is not same, any comparison of Fano factors could be misleading. We show, for example, that a change of intensity causes a false change of estimated Fano factor even if F remains fixed.

• This problem can be (mostly) eliminated by using Fano factor calculated in the operational time, when the values are compared for the greatest common operational time.

• If Fano factor is the same between different experiments, the operational Fano factor removes otherwise false change in the observed values. The bias of the differences of the estimated Fano factor is zero and MAE is (in studied cases) smaller than for the standard estimator.

• If the values of Fano factor change between the experiments, the situation is more complex. The operational Fano factor ensures better behavior of the estimates from the point of view of the influence of the intensity, however, it is still necessary to be careful about interpretation of the obtained values even for the operational Fano factor.

• When using operational Fano factor to compare various experiments, some data are wasted. This can be seen as a drawback of the method. It is, however, possible to use all the data so that we shift the shortened window and calculate the estimates of Fano factor multiple times and average them. We show that this procedure further decreases MAE of the estimated differences of Fano factor.

Data Availability Statement

Publicly available datasets were analyzed in this study. This data can be found at: http://dx.doi.org/10.6080/K0NC5Z4X.

Author Contributions

PL proposed the concept of the study. KR and PL wrote the manuscript. LK contributed the ideas and suggestions for improvement. KR performed the calculations, numerical simulations, and analysis. All authors contributed to manuscript revision and approved the submitted version.

Funding

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

Conflict of Interest

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.

References