## Abstract

Neuronal dynamics is intrinsically unstable, producing activity fluctuations that are essentially scale free. Here we study single cortical neurons of newborn rats *in vitro*, and show that while these scale-free fluctuations are independent of temporal input statistics, they can be entrained by input variation. Joint input–output statistics and spike train reproducibility in synaptically isolated cortical neurons were measured in response to various input regimes over extended timescales (many minutes). Response entrainment was found to be maximal when the input itself possesses natural-like, scale-free statistics. We conclude that preference for natural stimuli, often observed at the system level, exists already at the elementary, single neuron level.

## Introduction

Variability is a most prominent property of neural activity, both spontaneous and evoked. At the single neuron level, variability is observable in practically all aspects of evoked activity: irregularity of the spike train and trial-to-trial variability in spike counts, as well as irreproducibility of train structure evoked by identical input series (Faisal et al., 2008; Yarom and Hounsgaard, 2011). However, as demonstrated in several cases, response variability might be quenched by a variation introduced to the input itself (Churchland et al., 2010).

At the single neuron level, it was demonstrated that when stimulated with constant input, neuronal spike trains differ substantially between trials (Bryant and Segundo, 1976; Mainen and Sejnowski, 1995). In contrast, when stimulated with a fluctuating (filtered white noise) input, the reproducibility of the spike train is dramatically improved to the point of perfect repeatability, locking itself to (i.e., entrained by) input fluctuations, reliably encoding its structure. This key property was reproduced in a stochastic simulation of a Hodgkin–Huxley neuron, relating it to the properties of the underlying ion channels (Schneidman et al., 1998). While these measurements and simulations were limited to a timescale of seconds, it is known that when neuronal activity is observed over extended timescales, slower effects ensue (Marom, 2010). In a recent work (Gal et al., 2010) we have shown that, indeed, when presented with long (>1 h) sequences of pulse stimuli, single neuron response dynamics becomes intermittent and irregular, exhibiting scale-free fluctuations (e.g., with autocorrelation that lacks a characteristic scale). Given these slower modulatory processes, it is not obvious that the statistically unstructured random input series, which are capable of quenching response variation over limited timescales, will effectively entrain response variability over extended durations. The biophysical mechanism underlying the capacity of unstructured random input to entrain response variability relies on matching between timescales of input variations and timescales of the stochastic processes that generate the action potential. In contrast, when longer stochastic processes are allowed, they are left unmatched by the above unstructured input. It is therefore natural to hypothesize that to entrain neuronal response over extended timescales, the variations of the input series must match the scale-free temporal structure of intrinsic neuronal response dynamics. Indeed, at least at the system level, neuronal response variability is reduced under natural or natural-like sensory input (Aertsen and Johannesma, 1981; Baddeley et al., 1997; de Ruyter van Steveninck et al., 1997; Yu et al., 2005; Garcia-Lazaro et al., 2006, 2011). These natural-like signals are often characterized by long-range temporal correlations, and a general scale-free temporal structure (Voss and Clarke, 1975; De Coensel et al., 2003; Simoncelli, 2003).

In this study we directly measure the impacts of input temporal structure on response variability over extended timescales in isolated cultured cortical neurons. We show that while the response of neurons is temporally scale free, independently of input statistics, entrainment is maximal when the input itself has a matching, scale-free structure. We also perform analogous analyses to those of Mainen and Sejnowski (1995), quantifying the reproducibility of spike trains under different types of input. Here too, natural-like input minimizes the trial-to-trial variability of the spike train. We conclude that the rich and complex neuronal dynamics enable the neuron to match its dynamics to that of the natural environment, and that “tuning” to natural input statistics arises already at the atomic level of neural processing.

## Materials and Methods

##### Culture preparation.

Cortical tissues were obtained from newborn (<24 h) rats (Sprague Dawley) of either sex and dissociated following procedures described previously (Marom and Shahaf, 2002). The cells were plated directly onto substrate-integrated multi-electrode arrays (MEA) and developed for a time period of 2–3 weeks before their use. A total amount of ∼10^{6} cells was seeded on polyethyleneimine-coated MEAs. The preparations were kept in Minimal Essential Medium supplemented with heat-inactivated horse serum (5%), glutamine (0.5 mm), glucose (20 mm), and gentamycin (10 μg/ml), and maintained in an atmosphere of 37°C, 5% CO_{2}, and 95% air in an incubator as well as during the recording phases. An array of Ti/Au/TiN extracellular electrodes, 30 μm in diameter, and spaced either 500 μm or 200 μm from each other (MCS) were used. Synaptic transmission in the network was completely blocked by adding 20 μm amino-5-phosphonovaleric acid, 10 μm 6-cyano-7-nitroquinoxaline-2,3-dione, and 5 μm bicuculline-methiodide to the bathing solution. All experiments were performed in accordance with the regulations (and under the supervision) of the Technion–Israel Institute of Technology Animal Care Committee.

##### Measurements and stimulation.

A commercial amplifier (MEA-1060-inv-BC; MCS) with frequency limits of 150–3000 Hz and a gain of ×1024 was used. Rectangular 200 μs biphasic 600–800 mV voltage stimulation through extracellular electrodes was performed using a dedicated stimulus generator (STG4004; MCS). In the context of this study, no difference was observed in the behavior of neurons under current or voltage stimulation. Data was digitized to 16 bits using a USB-ME256 system (MCS). Each recorded channel was sampled at a frequency of 20 kHz. One hour after the addition of synaptic blockers, the stimulation electrode was selected as one evoking well isolated spikes with high signal-to-noise ratio, in as many recording electrodes as possible. From the selected recording electrodes, voltage traces of 15–20 ms poststimulus were collected. Spike detection was performed off-line by a manual threshold-based procedure. A 3 ms spike shape was extracted for each response for further noise cleaning and analysis. Stability of spike shape and activity dynamics criteria were applied to validate experimental stability, as described by (Gal et al., 2010). Random stimulation sequences were generated by modulating a constant stimulus interval sequence with a noise signal. This noise signal was either a white Gaussian noise or 1/*f* Gaussian noise generated by weighting the frequency components of the Gaussian white noise. In both cases a low cutoff was applied to have a minimum interval of 20 ms. The standard deviation (SD) of the noisy stimulation interval sequence was set such that its coefficient of variation will match the coefficient of variation of its response constant interval. For example, if a neuron responded to a 100 ms constant interval sequence with an average interspike interval (ISI) of 200 ms and SD of 40 ms, the SD of the noisy interval sequence was set to 20 ms.

##### Data analysis.

Analysis throughout this study was performed on either the spike time series, or on a smoothened firing rate time series, produced by filtering the spike train with a sliding rectangular window (Fig. 1*C*). Spectral analysis was performed on a binned time series, with a bin size of 1 s. The power spectrum density (PSD) was estimated using a modified periodogram. The Allan variance, which is commonly used to identify fractal point processes (Lowen and Teich, 1996, 2005) was calculated by binning the time series with different bin sizes *T*. For each binned time series, the Allan variance is defined as *Z _{T}* is the binned series, and <> mark the average over

*k*. Another measure, widely used to characterize the temporal statistics of long memory and nonstationary time series, is detrended fluctuation analysis (DFA; Peng et al., 1995), in which the fluctuations (in terms of mean square error) around a piecewise linear fit to the time series are quantified, for different segment durations. Since both the periodogram and the Allan variances can be regarded as power law for the purpose of this work, power law functions were fitted to the tail of these curves, and their exponents were used as descriptive statistics (Gal et al., 2010).

The similarity between two spike trains was assessed with two distance measures. (1) a correlation-based metric, defined as one minus the correlation between binned time series (*Z _{a}*,

*Z*):

_{b}*d*= 1 −

*corr*(

*Z*,

_{a}*Z*). This is a rate-based measure, which is insensitive to temporal features of the spike train below resolution dictated by the bin size. (2) The Victor—Purpura (VP) spike train metric (Victor and Purpura, 1997; Victor, 2005; Toups et al., 2011), which defines the distance between the two spike trains as the minimal cost of transforming one into the other. Briefly, a spike train is modified by a combination of three possible steps: inserting or deleting a spike, with a cost of +1, and moving a spike in time, with a cost of

_{b}*q*·

*dt*, where

*q*is a temporal resolution parameter. The value of

*q*sets the sensitivity of the metric to fine temporal features, and is set here to a default value of

*q*= 1

*s*. The dependence of the results on the value of

*q*is shown in Figure 4. The VP metric was calculated using MATLAB code downloaded from the Web site of Jonathan Victor.

## Results

### Unexplained response variability is minimized by scale-free input

In this paper we investigate the response properties of individual neurons, independent of synaptic and network effects. To that end, experiments are conducted on cultured cortical neurons, functionally isolated from their network by means of pharmacological block of both glutamatergic and GABAergic synapses (Gal et al., 2010; see Materials and Methods). Individual neurons are stimulated with sequences of short, identical extracellular electrical pulses. In response to a single pulse, a neuron either responds by emitting a spike, or fails to do so (Fig. 1*A*); neuronal responses are monitored by extracellular recording electrodes. As previously reported (Gal et al., 2010), when repeatedly stimulated over extended durations (minutes and more) with a low (∼1 Hz) stimulation rate, a neuron responds to each and every stimulation pulse (1:1 mode). As stimulation rate increases, the excitability of the neuron declines, and at a certain, neuron-specific critical stimulation rate, the 1:1 response mode breaks, and the neuron exhibits rich and complex response dynamics (Fig. 1*B*).

Here, we study the properties of response dynamics evoked by three statistically different regimes of stimulation series: (1) constant interval regime, (2) white noise regime in which the interval series is modulated by a Gaussian white noise process, and (3) scale-free regime in which the intervals are modulated by a 1/*f*^{β} process, with parameter β = 1, which is more representative of natural sensory input (Voss and Clarke, 1975; De Coensel et al., 2003; Simoncelli, 2003) and similar to activity properties of cortical neurons *in vivo* (Lowen and Teich, 1993, 1996; Teich et al., 1997; Lowen et al., 2001; Bhattacharya et al., 2005). The term scale-free is used here to designate a time series with an autocorrelation that decays slowly, usually as a power law without a typical scale. White noise on the other hand, which has a delta-function autocorrelation, is regarded as a zero scale signal. The choice of β = 1 is typical to a wide range of natural signals, in natural environments, and in biological systems. As will be shown, the results of this work are not sensitive to the exact value of β, consistent with previous studies (Garcia-Lazaro et al., 2006). Interval sequences in all stimulation regimes were normalized to have the same mean and SD (the latter is applicable only to the second and third regime). The mean stimulation rate (the reciprocal of the mean stimulation interval) was set to a high enough value to drive the neuron beyond its critical point, leading to response failures and intermittency (Gal et al., 2010). The interval SD was chosen to approximately match the intrinsic SD of the response to constant interval stimulation. Figure 1*D* shows examples of extracts from the three stimulation regimes, as well as their autocorrelation functions and the PSD.

Our analysis starts by stimulating neurons with long sequences (>1 h) of each of the stimulation regimes. Figure 2*A* shows an example of the response of the same neuron to the three regimes. It is immediately obvious that the three input regimes did not cause a significant difference in the statistical properties of the responses. This is formally shown in the plots of the PSD, Allan variance, DFA, rate histograms, and ISI histograms (Fig. 2*B–F*). Clearly, the macroscopic properties of the evoked spike trains are the same under all stimulation regimes: the neuron exhibits scale-free dynamics, characterized by power law statistics, in accordance with previously published analysis (Gal et al., 2010). To confirm the insignificance of the differences between the different statistical measures under the three regimes, a control experiment was performed, in which four neurons were stimulated and recorded, and each stimulation block was repeated 10 times. For each statistical analysis presented in Figure 2, *B–F*, the relevant parameter was estimated independently from the response to each repetition. The value ranges of these parameters from a single neuron are depicted in the insets of each part, and show that indeed the temporal statistics are the same under the three regimes. We emphasize that the claim made here is not that the specific model chosen for each statistic is the correct one (i.e., that the PSD is an exact power law), rather that these fits are good enough representative shapes, useful for comparing the population statistics. The above results imply that the various membrane and cellular processes underlying the stochastic response fluctuations are effectively insensitive to the statistical structure of the input regime.

1/*f*-type response statistics can be interpreted as a modulation of neuronal excitability by a cascade of oscillating processes at various timescales. An oscillator can be entrained (i.e., phase locked) by a driving stimulus at a frequency that matches the oscillator's natural frequency, and with a magnitude proportional to the oscillation amplitude. This suggests that the 1/*f* intrinsic fluctuations could be entrainable by a matching 1/*f* stimulation.

Figure 3 demonstrates this effect: the stimulation and response of a neuron are plotted in two timescales, for the three stimulation regimes. On the short timescale (Fig. 3*A*; 5 s bin size), the response in the white noise and scale-free regimes nicely follows the stimulation, while in the constant interval regime, the input has no bearing on timing of the output fluctuations. On a longer timescale (Fig. 3*B*; bin size of 100 s), the white noise input is practically flattened, becoming constant; as a result, it fails to entrain the response. In contrast, in the scale-free regime entrainment is evident throughout. Figure 3*C* shows a scatter plot of the input and output rates of the neuron, calculated with a 5 s bin size. As expected from the data of Figure 2, the input and output ranges of the white noise (red) and scale-free (green) regimes are practically identical. However, the correlation in the scale-free data is much higher, as expressed in the reduced variability around the linear trend line. This indicates that indeed the scale-free fluctuations in neuronal dynamics are best entrained by scale-free input. For 13 of the 17 recorded neurons, the input–output correlation coefficient significantly increases under scale-free input regime, compared with white noise regime. A Wilcoxon signed rank test yields *p* < 0.01 for the effect of the input regime on input/output correlation. Thus, under scale-free input, there is no change in the amount of variability in neuronal response, but more of it is explained by the input. Figure 3*D* shows a scaled correlation analysis: the correlation between input and output is calculated for responses in white noise and scale-free regimes on different timescales. The correlation in a given timescale *T* is calculated by smoothing the stimulation and response series with a rectangular window of width *T*, and subtracting from it a version of the series, smoothed with a window of size 25*T*. This effectively bandpasses the time series around *T*. While the correlation in the short timescale regime is similar for the white noise and scale-free regimes (the apparent increase for white noise is nontypical, a Wilcoxon signed rank test yields insignificant difference), for longer timescales it decreases in the case of white noise compared with scale-free regimes (*p* < 0.01, Wilcoxon signed rank test).

The above results show that there is a significant increase in correlation upon change from white noise (i.e., β = 0) to scale-free (β = 1) input. However, there is nothing unique about β = 1: other types of inputs might serve just as well. To assess the sensitivity of the increased input–output correlation to the value of β, neurons (*n* = 31) were exposed to nine stimulation blocks, each characterized by a different value of β, ranging from 0 to 2, and lasting 70 min. For each block, the correlation between input and response was computed and compared with the correlation for β = 0. The population statistics of this correlation ratio are depicted in Figure 3*E*, and show a clear concave shape, with a peak in the mid-range values, and a decrease toward lower and higher values of β. A similar peak is observed when plotting the distribution of preferred exponents (i.e., the exponents which results in the strongest correlation for each neuron; Fig. 3*F*). These peaks around β = 1, however, are wide, suggesting that the variability between neurons is considerable, and that inputs characterized by a wide range of exponents are equally effective in entraining neuronal responses.

The decrease of the correlation for input with large exponents, which are dominated by slow oscillations, is important for the understanding of this phenomena. Figure 3*G* shows extracts from the response to 1/*f* stimulation (left) and 1/*f*^{2} stimulation (right). It is easy to see how, for 1/*f*^{2}, the faster fluctuations remain unentrained, while slower oscillations are nicely locked by modulations in the input.

### Neuronal response repeatability is maximized by scale-free input

Another functionally relevant aspect of the entrainment described above concerns response repeatability. It is well established that despite the extensive variability of neuronal responses, spike train structure is repeatable when the input is fluctuating, both at the single neuron and sensory system levels (Mainen and Sejnowski, 1995; de Ruyter van Steveninck et al., 1997; Churchland et al., 2010). In the following set of experiments we ask whether a scale-free stimulation regime enhances repeatability in general, and over long timescales in particular. To this aim we stimulated neurons with 10 repetitions of the same input sequence under each of the three stimulation regimes: constant intervals, white noise, and scale free. Each sequence lasted 10 min, separated by a 10 min break. Figure 4*A* shows responses of one neuron to 10 identical sequences under each stimulation regime. It is immediately obvious that although some reproducibility exists under the white noise regime, the reproducibility of responses to scale-free sequences is much higher. This is a direct consequence of the previous section analysis: if responses are more correlated to the input, they will be more correlated between themselves under repetitions of the same input sequence. This observation is quantified using two kinds of spike train similarity measure: a rate-based measure (the correlation distance between spike histograms, calculated with 1 s bins; see Materials and Methods), and a time-based measure (Victor and Purpura, 1997; Victor, 2005; the VP distance, Toups et al., 2011). Figure 4*B* depicts the pairwise distance matrices of the responses of Figure 4*A*, calculated using these two metrics. Responses to white noise stimulation are significantly more reproducible than responses to constant input; this is in agreement with Mainen and Sejnowski (1995). But as expected from the results presented above (Fig. 3) the responses to the scale-free sequences are significantly more reproducible than those of white noise input. The purple lines of Figures 4, *C* and *D*, depict the mean and SD of the values in the distance matrices shown in Figure 4*B*. Of 14 recorded neurons, 11 showed significant improvement in reproducibility for scale-free sequences compared with constant interval and white noise, as quantified by the two metrics (there were no cases of disagreement). The gray lines depict the trends of mean distances calculated for all of the 14 recorded neurons. A Wilcoxon signed rank test for the mean pairwise distances shows an overall significant effect for the input regime on response repeatability (*p* < 0.01 for both white noise vs constant and for scale free vs white noise, for both metrics). The results are insensitive to the choice of the temporal parameter *q* of the spike train metric (Fig. 4*F*), which may point to the lack of characteristic scale in this phenomenon. The above enhanced reproducibility does not stem from a decreasing spike rate. Next to the preservation of the mean stimulation rate in all sequences, the spike rate itself does not decrease for scale-free input (Fig. 4*E*). For all of the 14 cells recorded, there was no significant decrease in firing rate under the scale-free regime compared with constant and white noise.

Since the scale-free input itself is structured, it is expected that even a “Bernoulli” neuron, which has a constant probability of response to a pulse stimulation, will have some reproducibility of its output spike train. Figure 5*A* shows a comparison between metric analyses on actual responses (gray), and on a surrogate Bernoulli neuron (purple) that responds with a constant probability (set to the mean response probability of the real neuron). As expected, the responses of the simulated neuron are indeed more reproducible for scale-free input, but not as consistent as the real neuron. It is also possible to construct a more detailed neuronal response curve, which takes into account the dependency of the response probability on the last interval between stimuli (an example for such a curve is given in Fig. 5*B*). Interestingly, the curves of the white noise and scale-free regimes substantially differ, pointing to a strong history dependence in response probability. As shown in Figure 5*A*, the resulting metric analyses (orange) behave more like the real neuron when compared with the Bernoulli neuron, yet cannot account for the entire effect. It is reasonable to believe that one might construct a response model that takes into account deeper history of the stimulation and response sequences to produce better fitting. It should be emphasized though, that while these neuronal response models do reproduce the repeatability effect to a significant extent, they do not reproduce the intrinsic scale-free fluctuations, and cannot be considered as successful explanatory models.

A scale-free input is characterized by the abundance of relatively long “breaks” in stimulation, or long periods with low stimulation rate, enabling recovery of internal processes from previous activations. As repeatedly shown over the past 15 years, the longer a neuron is exposed to repeated activations, longer recovery times are required (Toib et al., 1998; Ellerkmann et al., 2001; Fairhall et al., 2001; Lundstrom et al., 2008; Marom, 2009). Scale-free input statistics is inherently matched to such a mechanistic context: elongating a scale-free input series naturally gives rise to longer breaks, hence allowing for stabilization on every scale. An illustration of this property is provided by reanalyzing the data of Figure 4 over blocks of increasing lengths. The results are summarized in Figure 5, *C* and *D*, showing that the accumulation of variability (or divergence of response) with increasing block size (1–10 min range) is significantly slower than its accumulation in responses to other stimulation regimes.

## Discussion

In this paper we have shown how a natural-like, scale-free input entrains fluctuations of single neuron responses over extended timescales. We have demonstrated this property by comparing neuronal responses in three different stimulation regimes: constant interval, white noise, and scale-free. In the case of the scale-free regime, the correlation between the input and the response is significantly higher, and the repeatability of response is considerably enhanced. These characteristics are stable over long, practically unlimited durations. While the results do show a preference to mid-range values of β (around 1), there is nothing special about the exact value; what seems to be important is that the entrainment decreases when the slow-frequency component in the input becomes either too dominant (large β) or marginal (low β).

It has long been acknowledged that responsiveness of neural systems is optimized to the ranges of statistics found in natural inputs (Aertsen and Johannesma, 1981; Baddeley et al., 1997; de Ruyter van Steveninck et al., 1997; Yu et al., 2005; Garcia-Lazaro et al., 2006, 2011). Here we show that preference to natural statistics is not limited to large-scale neural systems; rather, it goes all the way down, to the atomic level of neural organization, namely the single, isolated neuron.

At the shorter timescales, the entrainment of neuronal fluctuations by white noise input is explained by the fast stochastic processes underlying the generation of action potentials (Schneidman et al., 1998). There, a variation in input allows for recovery of inactivation processes, unlike the case of a constant input that drives the neuron to operate around the limit of channel availability threshold, making it highly sensitive to stochastic events. It is reasonable to assume that a similar explanation would also be appropriate over extended timescales: a long break (or a period of low-rate stimulation) enables recovery of slow processes, in contrast to constant input (or shortly correlated input) that drives these processes to a highly stochastic operation point. Possible candidates for such processes might include slow inactivation properties of the ionic channels themselves (Toib et al., 1998; Ellerkmann et al., 2001; Marom, 2009; Soudry and Meir, 2010), or other cellular modulatory processes (e.g., protein phosphorylation, protein synthesis, and metabolic cycles).

From the more abstract, functional point of view, when the temporal structure of the input is relatively dull, it can only entrain a narrow range of cellular processes underlying neuronal dynamics. Under these conditions, a large fraction of the response variability is tagged “unexplained.” However, when the input is temporally rich, it matches the temporally complex structure of the intrinsic dynamics, and the former “unexplained” variability becomes information carrying. Thus, a neuron can be viewed as a collection of entangled information channels, distributed over a continuum of timescales. In this picture, information transfer is maximized when there is information to be transferred on any given scale.

## Footnotes

The research leading to these results has received funding from the European Union's Seventh Framework Program FP7 under grant agreement 269459 and was also supported by a grant of the Ministry of Science and Technology of the State of Israel and MATERA grant agreement 3-7878, and the ICRI-CI (Intel Collaborative Research Institute–Computational Intelligence). The work of A.G. was also partially supported by the Andrew and Erna Finci Viterbi Fellowship Program. We thank Danni Dagan, Daniel Eytan, and Avner Wallach for helpful discussions and insightful input, and Eleonora Lyakhov and Vladimir Lyakhov for invaluable technical assistance in conducting the experiments.

- Correspondence should be addressed to Asaf Gal, Network Biology Research Laboratories, Technion, 3200 Haifa, Israel. asaf.gal{at}mail.huji.ac.il