Broadband macroscopic cortical oscillations emerge from intrinsic neuronal response failures

Introduction

The most widespread cooperative activity of neurons within the cortex is spontaneous macroscopic oscillations (Silva et al., 1991; Gray, 1994; Contreras et al., 1997; Buzsáki and Draguhn, 2004; Chialvo, 2010), which range between sub- and hundred- Hertz (Başar et al., 2001; Brovelli et al., 2004; Buzsáki and Draguhn, 2004; Grillner et al., 2005; Giraud and Poeppel, 2012). The high cognitive functionalities of these oscillations are still controversial (Klimesch, 1996, 1999; Başar et al., 2001; Wiest and Nicolelis, 2003; Buzsaki, 2006; Kahana, 2006; Fries, 2009; Giraud and Poeppel, 2012; Thivierge et al., 2014) and are typically attributed to transitory binding activities among indirect macroscopic distant cortical regions (Gray, 1994; Başar et al., 2001; Buzsáki and Draguhn, 2004; Roxin et al., 2004; Fries, 2009). In addition, it was found that the theta rhythms (Klimesch, 1999; Buzsáki and Draguhn, 2004), oscillations in the range of 4-10 Hz, play a key role in the formation and retrieval of episodic and spatial memory (Hasselmo, 2005). This theta rhythm is usually accompanied by high frequency oscillations in the range of 30-80 Hz, known as gamma oscillations (Colgin and Moser, 2010). Gamma oscillations are also related to sensory stimulations and induce neuronal ensemble synchrony by generating a narrow window for effective excitation (Cardin et al., 2009).

There are several suggested mechanisms for the formation of such rhythms on the network level (Wang, 2010). Most of the proposed mechanisms are based on the existence of inhibitory synapses (Wilson and Cowan, 1972; Jirsa and Haken, 1996; Brunel and Wang, 2003), especially for high frequency oscillations (Brunel and Wang, 2003; Colgin and Moser, 2010; Wang, 2010). For illustration, assume that a fast excitation increases neural firing in an excitatory short-delayed feedback loop. Consequently, neuronal populations along the excitatory feedback loop will fire at high rates that will cause a slower response of the inhibitory neurons. As a result, the inhibitory neurons will depress the activity within the excitatory population. This will then depress the excitation of the inhibitory population. Finally, the depression of the inhibitory neurons allows a repeated fast excitation of the excitatory population.

In this work we show how extra-cellular potential oscillations, synchronized rhythmic firing of neurons, emerge in random networks without inhibitory synapses. Our findings are based on experimental observations of neuronal plasticity in the form of intrinsic neuronal response failures (NRFs; Vardi et al., 2015). Using simulations of large networks, based on single-neuron in-vitro experiments, we show that this type of neuronal plasticity leads to the coexistence of both theta and gamma oscillations. Results are supported by a quantitative approach based on a Langevin equation, which describes the network dynamics.

Materials and Methods

In-vitro Experiments

Animals

All procedures were conducted in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals and Bar-Ilan University Guidelines for the Use and Care of Laboratory Animals in Research and were approved and supervised by the Institutional Animal Care and Use Committee.

Culture Preparation

Cortical neurons were obtained from newborn rats (Sprague-Dawley) within 48 h after birth using mechanical and enzymatic procedures. The cortical tissue was digested enzymatically with 0.05% trypsin solution in phosphate-buffered saline (Dulbecco's PBS) free of calcium and magnesium, and supplemented with 20 mM glucose, at 37°C. Enzyme treatment was terminated using heat-inactivated horse serum, and cells were then mechanically dissociated. The neurons were plated directly onto substrate-integrated multi-electrode arrays (MEAs) and allowed to develop functionally and structurally mature networks over a time period of 2–3 weeks in-vitro, prior to the experiments. Variability in the number of cultured days in this range had no effect on the observed results. The number of plated neurons in a typical network was in the order of 1,300,000, covering an area of about 380 mm². The preparations were bathed in minimal essential medium (MEM-Earle, Earle's Salt Base without L-Glutamine) supplemented with heat-inactivated horse serum (5%), glutamine (0.5 mM), glucose (20 mM), and gentamicin (10 g/ml), and maintained in an atmosphere of 37°C, 5% CO₂ and 95% air in an incubator as well as during the electrophysiological measurements.

Synaptic Blockers

All experiments were conducted on cultured cortical neurons that were functionally isolated from their network by a pharmacological block of glutamatergic and GABAergic synapses. For each culture 20 μl of a cocktail of synaptic blockers was used, consisting of 10 μM CNQX (6-cyano-7-nitroquinoxaline-2,3-dione), 80 μM APV (amino-5-phosphonovaleric acid), and 5 μM bicuculline. This cocktail did not block the spontaneous network activity completely, but rather made it sparse. At least 1 h was allowed for stabilization of the effect.

Stimulation and Recording

An array of 60 Ti/Au/TiN extracellular electrodes, 30 μm in diameter, and spaced 500 μm from each other (Multi-Channel Systems, Reutlingen, Germany) was used. The insulation layer (silicon nitride) was pre-treated with polyethyleneimine (0.01% in 0.1 M Borate buffer solution). A commercial setup (MEA2100-2 × 60-headstage, MEA2100-interface board, MCS, Reutlingen, Germany) for recording and analyzing data from two 60-electrode MEAs was used, with integrated data acquisition from 120 MEA electrodes and 8 additional analog channels, integrated filter amplifier and three-channel current or voltage stimulus generator (for each 60 electrode array). Mono-phasic square voltage pulses typically in the range of [−800, −500] mV and [60, 400] μs were applied through extracellular electrodes. Each channel was sampled at a frequency of 50 k samples/s, thus the changes in the neuronal response latency (NRL) were measured at a resolution of 20 μs.

Cell Selection

Each node was represented by a stimulation source (source electrode) and a target for the stimulation—the recording electrode (target electrode). These electrodes (source and target) were selected as the ones that evoked well-isolated, well-formed spikes and reliable response with a high signal-to-noise ratio. This examination was done with a stimulus intensity of −800 mV with duration of 200 μs using 30 repetitions at a rate of 5 Hz, followed by 1200 repetitions at a rate of 10 Hz.

In-vitro Experiment with Feedback Loops and Neural Circuits

The activity of all source and target electrodes was collected and action potentials were detected on-line by threshold crossing, and entailed stimuli were delivered in accordance with the circuit's connectivity, as described below. A successful response was defined as a spike occurring within a typical time window of 2–10 ms following the beginning of an electrical stimulation.

In Figures 2A,B, after every spike detection two supra-threshold extracellular stimulations were given to the same neuron, after 600 ms and 630 ms. In case that the timings of the stimulations overlap, only one stimulation is given.

In Figures 2C,D, after every spike detection supra-threshold extracellular stimulations were given to its connected neurons. For example, if a spike was detected at the left (green) neuron (Figure 2C), a supra-threshold extracellular stimulation will be given to the middle (brown) neuron after 330 ms.

Data Analysis

Analyses were performed in a Matlab environment (MathWorks, Natwick, MA, USA). The reported results were confirmed based on at least eight experiments each, using different sets of neurons and several tissue cultures.

The temporal firing frequency, around stimulation no. i, of the neuron in Figure 1 was calculated using the following procedure

$\begin{array}{l}firingfrequency\left(i\right)=stimulationfrequency·\\ {\displaystyle \sum _{\text{m=max}(0,i-125)}^{i+125}\frac{Is\_Spiked(m)}{i+125-\max (0,i-125)}}\end{array}$

where Is_Spiked(m) = 1 if the neuron responded to stimulation no. m, otherwise Is_Spiked(m) = 0.

FIGURE 1

Figure 1. Neuronal plasticity—in-vitro experiments. Upper panel: The NRL of a neuron stimulated at 10, 12, and 15 Hz (light blue, green, and purple dots, respectively). Response failures are denoted by NRL < 2 ms. Lower panel: The firing frequency calculated from the averaged ISI using sliding windows of 250 stimulations, or the maximal available one for Stim < 250.

Simulations

Simulations (similar to Vardi et al., 2015) consist of a network of N leaky integrate and fire neurons

$\frac{d{V}_i}{dt}=-\frac{V_i}{\tau }+{\displaystyle \sum _{j=1}^N{J}_{ji}{\displaystyle \sum _{t^{\prime }}\delta \left(t-t\prime -D_{ji}\right)}}$

where i∈[1, N], τ = 20 ms, J_ji and D_ji are the connection's strength and delay from neuron j to i, respectively. The summation over t′ sums all firing times of neuron j, the integration time step is 0.05 ms, and the threshold is 1. For the n^th threshold crossing of a neuron, its probability for a response is [Σ(τ_n−m/τ_c)exp(-αm)]/[Σexp(-αm)], where τ_n is the time gap between the n^th and (n–1)^th threshold crossings, τ_c = 1/f_c, α = 1.4 and the sum is over m≥0. A refractory period of 2 ms is imposed after an evoked spike, for response failures the voltage is set to 0.2. Results were found to be insensitive to initial conditions.

Various Forms of p(s|τ) Lead to the Same <ISI>

The probability for a response, given the last inter-stimulation-interval, p(s|τ), should lead to < ISI> = τ_c (Figures 1, 3). One can show that any p(s|τ) satisfying

$\frac{{\displaystyle {\int }_0^{\infty }p\left(s|\tau \right)p\left(\tau \right)d\tau }}{{\displaystyle {\int }_0^{\infty }\tau p\left(\tau \right)d\tau }}=\frac{1}{{\tau }_c}$

where p(τ) is the probability density of an inter-stimulation-interval equals to τ, leads to < ISI> = τ_c. The numerator on the left hand side stands for the average probability for a successful response, and the denominator stands for the average inter-stimulation-interval. This ratio is equivalent to the firing rate, hence equals to 1/τ_c. For instance p(s|τ) = τ/τ_c, this theoretical curve fits all p(τ) (Figure 3D). In the activity of some random networks a good approximation for p(τ) is

$p\left(\tau \right)=\frac{2}{{\tau }_c}exp\left(-\frac{2\tau }{{\tau }_c}\right).$

For this p(τ) some of the p(s|τ) solutions, which lead to < ISI> = τ_c, are p(s|τ) = 0.5, p(s|τ) = τ/τ_c and

$p\left(s|\tau \right)=1-exp\left(-\frac{2\tau }{{\tau }_c}\right),$

which is similar to Figure 3D.

Fourier Analysis of the Rate

To perform a Fourier analysis on the activity of the network we define the rate vector:

$R\left(time=i·w\right)=\frac{1}{Nw}{\displaystyle {\int }_{(i-0.5)w}^{(i+0.5)w}{\displaystyle \sum _{t_{spike}}\delta \left(t\prime -t_{spike}\right)}dt\prime }$

where i is an integer, w is a predefined time window and the sum is over all spike times of all N neurons.

Next, a discrete Fourier transform is preformed and the function of resulting amplitudes is normalized and smoothed using a sliding window of 1 Hz.

Results

Neuronal Plasticity: Neuronal Response Failures

We start with the quantification of the NRL (Vardi et al., 2013a, 2015; Marmari et al., 2014), measured as the time-lag between a stimulation and its corresponding evoked spike (Figure 1). It was recently shown (Vardi et al., 2015) that when a neuron is repeatedly stimulated, its NRL stretches gradually (Figure 1, upper panel), and when the stimulation frequency is high enough stochastic NRFs emerge. Specifically, each neuron is characterized by a critical frequency, f_c, typically ranging among neurons between 2 and 30 Hz. Stimulation frequencies above f_c result in NRFs, whereas for supra-threshold stimulations below f_c a response is assured. The probability of the NRFs is such that the neuron functions similar to a low pass filter, saturating its firing rate (Figure 1, lower panel). Quantitatively, for a stimulation frequency f, the NRF probability is 1–f_c/f, i.e., the firing frequency is saturated at f_c. Thus, changing the stimulation frequency will change the probability for NRFs, while the firing frequency remains bounded, f_c.

This observation is demonstrated using a cultured cortical neuron, functionally separated from its network by synaptic blockers, with above-threshold stimulations (see Section In-vitro Experiments). We examine a neuron with periodic stimulation trials of 10, 12, and 15 Hz, and NRFs appear after a short transient where the neuron exhibits an increase of its NRL (Figure 1, upper panel). Examining the temporal firing rate of the neuron (Figure 1, lower panel), it is noticeable that the firing frequency of the neuron is saturated at f_c = 5 Hz, independent of the stimulation frequency.

The effect of NRFs is first examined on small neuronal circuits using the following experiment: We stimulate the neuron ones and impose on the neuron two self-feedback delay loops, e.g., 600 and 630 ms (Figure 2A). The neuron is stimulated 600 and 630 ms after each evoked spike (see Section In-vitro Experiment with Feedback Loops and Neural Circuits). In the case of vanishing probability for NRF, the neuron should fire every 30 ms (Kanter et al., 2011; Vardi et al., 2012b), i.e., 33.3 Hz. Since f_c = 5 Hz, some NRFs appear (Figure 2B) forming bunched firing, separated by ~600 ms, whereas the intra-bunch time-lags, 30 ms, originated from the difference between the two feedback delay loops (Figure 2B). The emergence of such firing bunches indicate some dynamical changes of the NRF's probability, where the neuron adapts its failure probability as a result of its recent stimulation history.

FIGURE 2

Figure 2. Firing bunches stem from neuronal response failures—in-vitro experiments. (A) Schematic of the neuron in Figure 1 characterized by f_c = 5 Hz, with 600 and 630 ms self-feedback loops. (B) A 3.5 s snapshot of the experimental results of (A) where stimulations (blue lines) and their corresponding evoked spikes (blue dashed lines) were recorded after an offset of t_o = 21 s and a preparation at 5 Hz stimulation frequency over 300 s to settle the neuron at the intermittent phase. (C) Schematic of a circuit consisting of three different neurons with f_c = 6.8 (green), 7.2 (brown), and 4.0 (purple) Hz and 600 and 630 ms delay loops, similar to (A), [different neurons than the one in (A)]. (D) A 5 s snapshot of the experimental results of (C) where 10 initial stimulations were given to the central neuron (brown) at 4 Hz and t_o = 25 s. Stimulations given to the colored neurons (brown/green/purple lines, respectively) and their corresponding evoked spikes (brown/green/purple dashed lines, respectively). A zoom-in of the gray area is presented (bottom).

A more biologically realistic scenario is a neuronal circuit consisting of three artificially connected neurons (Vardi et al., 2012b), forming the same delay loops (Figure 2C), 600 and 630 ms. In addition, each neuron is identified by different f_c and only the central neuron receives the initial stimulation. For the central neuron, which is characterized by f_c = 7.2 Hz, NRFs occur, since the frequency of its driven neurons is greater than its critical frequency, 6.8+4>7.2 Hz, and similarly the outer neurons have response failures (7.2 > 6.8, 7.2 > 4) (Figure 2D). Besides the formation of firing bunches for each neuron, their firing are correlated at zero- or shifted- time-lags (Figure 2D). The question whether the repeated bunches in such small neuronal circuits shed light on macroscopic cortical oscillations requires to sail toward large scale simulations.

Short-term Memory of Neuronal Plasticity

The observed firing bunches indicate a form of short-term memory of neuronal plasticity, where the NRF probability is mainly a function of the few preceding stimulations. Our next goal is to experimentally quantify this neuronal plasticity and then examine its implementation on the dynamics of large neural networks using large scale simulations.

We first assume a simplified network where each neuron has two above-threshold inputs, two outputs and the same f_c (Figure 3A), hence the statistics of the inter-stimulation-intervals of each neuron is expected to approximately follow an exponential distribution with 2f_c rate, Exp(2f_c) (Figure 3B, upper panel). To quantify the statistics of the NRFs, a long stimulation trail obeying Exp(2f_c), under τ_c/8 time resolution, was given to a cultured neuron with f_c ~ 5 Hz (Figure 3B). Next, the conditional probabilities for a successful response (an evoked spike), p(s|τ_i,τ_i−1), were estimated for events where the current inter-stimulation-interval equals τ_i and the previous one equals τ_i−1 (Figures 3B,C). It is evident that p(s|τ_i,τ_i−1) is primarily a function of τ_i, i.e., the probability for a successful response is dictated mainly by the current inter-stimulation-interval, τ_i. Hence, the NRFs can be dynamically approximated using p(s|τ) (Figure 3D), indicating that the NRF probability might be fairly estimated based solely on the last inter-stimulation-interval, τ.

FIGURE 3

Figure 3. The short-term memory of the stochastic neuronal response failures—in-vitro experiments and simulations. (A) Schematic of a prototypical examined excitatory network, where each neuron has two pre- and two post- synaptic connections and the same f_c. Inset: The time-lags between neuronal stimulations is expected to approximately follow an exponential distribution, Exp(2f_c). (B) Upper panel: The stimulation scheme where a neuron is stimulated under the resolution of τ_c/8, such that the discrete differences between two consecutive stimulations, τ, follows Exp(2f_c). Middle panel: Experimental NRL of a cultured neuron with f_c ~ 5 Hz under a long trial of stimulations following Exp(2f_c), response failures are denoted at NRL = 3.5 ms. Lower panel: Zoom-in of the middle panel (gray area) and schematic of the conditional probability p(s|τ_i,τ_i−1), measuring the probability of a successful response, spike, given that the current inter-stimulation-interval equals τ_i and the previous one equals τ_i−1. (C) The probabilities p(s|τ_i,τ_i−1) obtained from the experimental data in (B) for time >3500 (middle panel). Each colored line presents p(s|τ_i,τ_i−1) for a given τ_i−1 in τ_c/8 time units (legend). (D) The experimentally measured p(s|τ_i) (black), measuring the probability of a successful response, spike, given that the current inter-stimulation-interval equals τ_i, and the theoretically predicted one (green) using the simplified assumption, p(s|τ_i) = τ_i/τ_c for τ_i < τ_c. For both curves, the average ISI ~ τ_c is preserved. (E) The probability density function of inter-stimulation-intervals, τ, for all neurons of Figure 4A (blue). (F) Typical Fourier amplitude of spike timings of a randomly chosen neuron, taken from Figure 4A.

Neuronal Oscillations on the Network Level

The experimental estimation of p(s|τ) is utilized to simulate a large scale network (Figure 3A) and is exemplified for 2000 neurons where delays between connected neurons, D, are randomly selected from a uniform distribution U(10, 15) ms. The simulation is initialized with timings of evoked spikes for a subset of the neurons, however, besides the transient time results were found to be insensitive to the initial conditions. The response failure is then randomly selected following the experimentally measured p(s|τ), independently for each neuron and stimulation (Figure 3D). Indeed, the assumption Exp(2fc) (Figure 3B, upper panel) was confirmed (Figure 3E), the statistics of the inter-stimulation-intervals of each neuron approximately follows an exponential distribution with 2f_c rate. The raster plot of the network firing as well as the time-dependent firing rate (Figures 4A,B) clearly indicate cooperative oscillation which can be quantified using the Fourier analysis to f_osc ~ 3 Hz (Figure 4C, see Section Simulations), and are absent in the Fourier analysis of the firing of each individual neuron (Figure 3F). Another broadened Fourier peak is centered at f_γ ~ 80 Hz, gamma oscillations (Brunel and Wang, 2003; Cunningham et al., 2004; Fries, 2009; Minlebaev et al., 2011; Dugladze et al., 2012), which is attributed to the average delay, < D> = 12.5 ms. It reflects the average firing frequency of each neuron where NRFs vanish and all delays are equal to < D>, as GCD = < D> for delay loops of such a random network (Kanter et al., 2011; Vardi et al., 2012a). Similar cooperative oscillations were obtained using a counterpart simulation for the same network (Figures 4D,F) while using the theoretical form p(s|τ) (Figure 3D). Although the form of p(s|τ) varies among neurons as well as between the theoretical and the experimental forms (Figure 3D), the cooperative oscillations are found quantitatively to be only slightly affected by its exact form. The robustness of f_osc was also confirmed in simulations for more realistic scenarios where f_c significantly varies among neurons as well as their input and output connectivity distributions. A more biological realization is exemplified in Figures 4G–I. The number of connections per neuron is much greater than 2, i.e., more than 50 pre- and 50 post- synaptic connections, where most of them are sub-threshold and on the average 1.5 of pre- and post-synaptic are above-threshold. Specifically, each sub-threshold connection produces an excitatory postsynaptic potential, EPSP, which is equal to 0.03 of the threshold. It is apparent that these additional connections do not qualitatively change the oscillations. In addition, f_γ was found to be robust to a wider distribution of delays and followed its center (Figure 4I). Note that without these intrinsic NRFs, i.e., p(s|τ) = 1, the firing of each neuron is only bounded by the refractory period which is in the order of several milliseconds. In this limiting case, the abovementioned theta and gamma oscillations disappear, as was confirmed in simulations (not shown).

FIGURE 4

Figure 4. Cooperative cortical oscillations on a network level. (A) Raster plot of the evoked spikes (blue dots) obtained in the simulation of a network of 2000 neurons with f_c = 5.7 Hz. Each neuron has two randomly selected pre- and post- synaptic connections, and the simulation is based on the experimentally obtained p(s|τ), (Figure 3D). Delays are randomly selected from U(10, 15) ms. The contrast of the raster was enhanced using a dilution of a constant amount of randomly chosen points in each sliding window of 23 ms, with a step of 0.23 ms. The average dilution is ~60% of the points. (B) The average firing rate per neuron as a function of time, calculated for windows of 1 ms. (C) The normalized Fourier amplitude, using a sliding window of 1 Hz, of the entire firing of all neurons over a time slot of 30 s, indicating f_osc ~ 3.6 Hz and f_γ ~ 80 Hz. Inset: The normalized Fourier amplitudes in the range [0,30] Hz obtained from R(m), Equation (1), D = 12.5 ms (red) and from the simulation, [(B), blue]. (D–F) Similar to (A–C) where each neuron has on average two pre- and post- randomly chosen synaptic connections and utilizing the theoretical p(s|τ), (Figure 3D). f_osc ~ 4.0 Hz and f_γ ~ 80 Hz at (F). (G–I) Simulation of a network of 2000 neurons where each neuron has on the average 1.5 pre- and 1.5 post- synaptic above-threshold connections, and 50 pre- and 50 post- synaptic sub-threshold connections with a strength of 0.03, relative to a threshold of 1. p(s|τ_i) is generalized to an exponential decay function of the neuronal stimulation history, Σ(τ_i−m/τ_c)exp(-αm)/[Σexp(–αm)], α = 1.4 and the sum is over stimulation history, m≥0. Delays are randomly selected from U(12.5, 20) ms and f_c from U(3, 10) Hz. f_osc ~ 2.6 Hz and f_γ ~ 60 Hz at (I), and the inset is similar to (C) and (F), but with K = 1.5, f_c = 6.5 Hz and D = 16.25 ms.

Cortical Oscillations vs. Statistical Properties of the Network

The origin of the fast oscillations, f_γ, (Figures 4C,F,I) is evident, however, the mechanism underling the slow cooperative oscillations (Wu et al., 1999; Sanchez-Vives and McCormick, 2000; Bollimunta et al., 2008, 2011; Nir et al., 2008; Crunelli and Hughes, 2010), f_osc, is still unclear. To explore this mechanism we identify the following two characteristic distances on the network. The first distance, Path, is the average minimal path between pairs of nodes, and the second distance, Loop, is the average over the minimal feedback loop of each node, both counted by the number of nodes along the route (Figure 5A). Numerical estimations based on various network topologies indicate that the average values of these two distances as well as their distributions are almost identical (Figure 5B) and their scaling decrease as 1/ln(< K>), where < K> is the average neuronal input connectivity (Figure 5C). These identities and scaling (Figures 5B,C), are also supported by the following theoretical argument. Assume a random network consisting of N neurons and an average connectivity < K>. The quantity Q(m) denotes the number of new connected neurons to a seed neuron at a distance of m neurons, hence Q is proportional to the probability (green) in Figure 5B. Start at an arbitrary neuron, Q(0) = 1, this neuron is connected (pre-synaptic) to Q(1) neurons. Using a recursive formula one can approximate

$\text{Q}\left(\text{i}\right)=\text{N}\left(\text{1}-\text{exp}\left(-<\text{K}>\frac{\text{Q}\left(\text{i}-\text{1}\right)}{\text{N}}\right)\right)\left(1-{\displaystyle \sum _{\text{m}=\text{1}}^{\text{i}-\text{1}}\frac{\text{Q}\left(\text{m}\right)}{\text{N}}}\right)$

where N(1−exp(− < K>Q(i–1)/N)) stands for the average number of neurons at a distance i, which are connected from new neurons at distance i–1. The rightmost term is the probability that the neuron at distance i is a new neuron which was not counted at shorter distances, m < i. This recursive relation is solved numerically and the normalized Q(i) are presented in Figure 5B. The above analysis is valid for Loops and Paths, hence their statistics are identical, in agreement with the sampling of these quantities in Figure 5B.

FIGURE 5

Figure 5. Scaling properties of f_osc and f_γ. (A) Schematic of an excitatory network where a self-feedback loop (light-red line) and the minimal self-feedback loop (red line) for a given neuron (filled red circle) are denoted. Similarly, a path (light-blue line) between two neurons (filled blue circles) and the minimal path (blue line) are denoted. (B) The distribution and its average (vertical lines) for the minimal path (Path) and for the minimal loop (Loop) obtained in simulations for networks as in (Figure 3A) with N = 4000, error bars are comparable with the circles. The analytical estimation is shown in green. (C) The scaling of the averaged quantities in (B) as a function of the average connectivity, < K>. (D) Simulation results indicate f_osc ∝ ln < K>, where N = 4000, f_c is randomly chosen for each neuron from U(5, 15) Hz and delays are randomly chosen from U(10, 15) ms. The probability for a connection between two neurons is < K>/N. Error bars indicate the standard deviation. (E) Simulation results indicate f_osc ∝ln(f_c), for networks as in (D), with < K> = 2, but f_c is the same for all neurons. (F) Simulation results indicate f_γ ∝ 1/ < D>, for networks as in (D), with < K> = 2, but delays are randomly chosen from U(< D> −2.5, < D>+ 2.5) ms.

The importance of this distance in the formation of f_osc can be understood according to the following argument. Assume a random subgroup of firing neurons activates another random group of neurons and vice versa. The minimal delay between pairs of neurons belonging to the two subgroups is Path· < D>. Consequently, the oscillations are expected to scale with Path· < D>, and indeed results indicate that f_osc∝ ln(< K>) ∝ (Path)⁻¹ (Figures 5C,D). The minimal path is the most reliable one with respect to the NRFs, however, the effect of longer paths is not negligible as they might maintain the activity of the Path (Figure 2), especially as their entropy is higher. In addition, the NRFs are responsible to limit the firing of all the network simultaneously (Figures 4B,E,H). Assume one neuron fires and activates < K> neurons after < D> ms, hence after m < D> ms, < K>^m neurons fire. This exponential firing growth is bound by m ≈ Path, as self-feedback loops (Figures 2A,B) significantly lead to NRFs and to a fast decrease in the firing rate. In addition, for a given network topology, f_osc is found to scale with ln(f_c) (Figure 5E), and to be robust for networks composed of neurons with different f_c (Figures 4G,H,I). These predictions might be realized in further experiments by controlling the network topology either by the neuronal concentration or by pharmacological manipulations.

Analytical Description of the Network Oscillations

An analytical description of f_osc is also possible, and to simplify the presentation the method is briefly described for homogeneous networks only. Each node has the same fixed input and output connectivity, K, and all delays are equal to D ms, hence neurons can fire only at i·D ms, where i is an integer. The fraction, R(m), of neurons that fire at step m is given by

$\begin{array}{cc}R\left(m\right)=\chi \left(m\right)·K·R\left(m-1\right)+\xi \left(m\right)& \left(1\right)\end{array}$

where x(m) stands for the time-dependent white noise representing the stochastic nature of the response failures. The function c(m) represents the susceptibility of the network, i.e., the fraction of neurons that fires if all neurons are stimulated at step m, and is explicitly given by

$\begin{array}{cc}\chi \left(m\right)={\displaystyle \sum _{i=0}^{\infty }p(s|i)·K·R(m-\left(i+1\right)){\displaystyle \prod _{n=2}^i\left(1-K·R(m-n)\right)}}& \left(2\right)\end{array}$

The first term, p(s|i), stands for the probability for an evoked spike when the previous stimulation was given before i steps (Figure 3D). The term K·R(m–(i + 1)) pinpoints a neuron stimulated before i steps, and the product, the rightmost term, indicates that the neuron was not stimulated since step (m–i). Equation (2) indicates that χ(m) is a function of the variable R(l) only, with l<m. Hence, after the insertion of Equation (2), into Equation (1), one finds a recursion relation for R(m) which can be solved numerically given the initial conditions. The dynamical solution of this recursion relation revealed oscillations which were found to fit fairly good with those observed in large-scale simulations (Figures 4C,F,I). The equations imply that the network has some memory of its previous activity, which dictates the responsiveness of the entire network. This analytical description can be generalized to advanced structured networks, including random connectivity, distribution for the delays as well as to include variations among neuronal critical frequencies, f_c.

Discussion

We have demonstrated that intrinsic NRFs drive a neural network activity toward oscillations, where high frequency oscillations, gamma, and low frequency oscillations, delta and theta, coexist. The high frequency oscillations correspond to the average delay between connected neurons in the network, while low frequency oscillations are governed by statistical properties of the network, e.g., the average number of connections per neuron and the average critical frequency of neurons. The coexistence of high and low frequency oscillations was confirmed in a new type of simulations, based on a single neuron in-vitro experiments, to evaluate the firing activity of complex networks. Results were also supported by an analytical description of the stochastic dynamics of the network.

Preliminary results in-vivo support our findings. The NRL increases by several milliseconds under periodic stimulations and terminates at an intermittent phase (Vardi et al., 2015). This phase is characterized by fluctuations around a constant NRL and accompanied by NRFs. Results indicate that f_c can be below 10 Hz and vary among neurons. However, quantitative measurements of f_c and the statistics of the NRFs require long stimulation trials, i.e., many thousands of high frequency periodic stimulations, which are currently beyond our experimental capabilities.

The average firing rate of neurons in the network is low, e.g., ~3.6, ~3.9, and ~2.6 Hz in Figures 4B,E,H, respectively. These network low firing rates are lower than the neuronal critical frequency in Figures 4B,E, f_c = 5.7 Hz, and < f_c> = 6.5 Hz, in Figure 4H. Surprisingly, the neuronal critical firing frequency is not saturated even when the network is completely excitatory. A biological mechanism that suppresses the firing frequency of a single neuron below f_c is aperiodic time-lags between stimulations (Vardi et al., 2015). For illustration, assume a slow mode of alternation between stimulation frequencies of 2f_c (0.5τ_c) and 2f_c/3 (1.5τ_c), such that < τ> = τ_c. For the high and low frequency mode, the expected probability for a NRF is 0.5 and 0, respectively. Consequently, < ISI> = 0.5(1.5τ_c + τ_c) = 1.25τ_c, corresponding to a lower firing rate, 0.8f_C. In addition, the firing rates are considerably lower than f_γ, indicating that high frequency network oscillations consist of temporarily synchronized sub-groups of neurons. Indeed, the Fourier spectrum of a single neuron does not exhibit any dominant peaks (Figure 3F).

Although the formation of broadband network oscillations is usually attributed to the existence of inhibition, it is evident that another possible mechanism is intrinsic NRFs that dynamically drives neural networks to generate coherent oscillations with low averaged firing rates (Vardi et al., 2015). These observations raise the question of which functionalities demand synaptic inhibition. It was shown that inhibition slightly suppresses the network firing frequency even further (Vardi et al., 2015) and it also might change the amplitude of the oscillations. An additional possible hypothesis is that the role of inhibitory connections is to allow some neuronal computations which are based on conditional temporal formation of neuronal firing patterns. This type of functionality is an exclusive property of inhibitory synapses which probabilistically block an evoked spike of its driven nodes in a given time window (Vardi et al., 2013b; Goldental et al., 2014). In addition, the coexistence of the network oscillations with neuronal inhibition is intriguing, and especially the questions whether inhibition induces more modes of oscillations, sharpens the existing ones, or suppresses the oscillatory behavior and stabilizes the network activity. Preliminary results of simulations indicate that inhibition might suppress the amplitude of the oscillations in the low frequency range and sharpen the oscillations in the gamma range. However, results might be sensitive to the selected parameters.

The interplay between the presented spontaneous cortical oscillations and external stimulations given to the network is another intriguing question. Specifically, it is interesting to examine the coexistence and the interplay between the spontaneous oscillation frequencies determined by the network topology and the frequencies of the induced external stimulations. The understanding of these dynamics will shed light on the emerging cortical oscillations among coupled networks characterized by different statistical properties.

Author Contribution

RV and SS prepared and performed the experiments; RV, SS, and AG analyzed the data; PS designed the experimental real-time interface; AG performed the simulations and developed the theoretical framework with help of PS.; IK, RV, and AG wrote the manuscript.

Conflict of Interest Statement

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.

Acknowledgments

This research was supported by the Ministry of Science and Technology, Israel. The authors declare no competing financial interests.

References

Başar, E., Baar-Eroglu, C., Karaka, S., and Schürmann, M. (2001). Gamma, alpha, delta, and theta oscillations govern cognitive processes. Int. J. Psychophysiol. 39, 241–248. doi: 10.1016/S0167-8760(00)00145-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Bollimunta, A., Chen, Y., Schroeder, C. E., and Ding, M. (2008). Neuronal mechanisms of cortical alpha oscillations in awake-behaving macaques. J. Neurosci. 28, 9976–9988. doi: 10.1523/JNEUROSCI.2699-08.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Bollimunta, A., Mo, J., Schroeder, C. E., and Ding, M. (2011). Neuronal mechanisms and attentional modulation of corticothalamic alpha oscillations. J. Neurosci. 31, 4935–4943. doi: 10.1523/JNEUROSCI.5580-10.2011

PubMed Abstract | CrossRef Full Text | Google Scholar

Brovelli, A., Ding, M., Ledberg, A., Chen, Y., Nakamura, R., and Bressler, S. L. (2004). Beta oscillations in a large-scale sensorimotor cortical network: directional influences revealed by Granger causality. Proc. Natl. Acad. Sci. U.S.A. 101, 9849–9854. doi: 10.1073/pnas.0308538101

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunel, N., and Wang, X.-J. (2003). What determines the frequency of fast network oscillations with irregular neural discharges? I. Synaptic dynamics and excitation-inhibition balance. J. Neurophysiol. 90, 415–430. doi: 10.1152/jn.01095.2002

PubMed Abstract | CrossRef Full Text | Google Scholar

Buzsaki, G. (2006). Rhythms of the Brain. New York, NY: Oxford University Press. doi: 10.1093/acprof:oso/9780195301069.001.0001

CrossRef Full Text | Google Scholar

Buzsáki, G., and Draguhn, A. (2004). Neuronal oscillations in cortical networks. Science 304, 1926–1929. doi: 10.1126/science.1099745

PubMed Abstract | CrossRef Full Text | Google Scholar

Cardin, J. A., Carlén, M., Meletis, K., Knoblich, U., Zhang, F., Deisseroth, K., et al. (2009). Driving fast-spiking cells induces gamma rhythm and controls sensory responses. Nature 459, 663–667. doi: 10.1038/nature08002

PubMed Abstract | CrossRef Full Text | Google Scholar

Chialvo, D. R. (2010). Emergent complex neural dynamics. Nat. Phys. 6, 744–750. doi: 10.1038/nphys1803

CrossRef Full Text | Google Scholar

Colgin, L. L., and Moser, E. I. (2010). Gamma oscillations in the hippocampus. Physiology 25, 319–329. doi: 10.1152/physiol.00021.2010

PubMed Abstract | CrossRef Full Text | Google Scholar

Contreras, D., Destexhe, A., Sejnowski, T. J., and Steriade, M. (1997). Spatiotemporal patterns of spindle oscillations in cortex and thalamus. J. Neurosci. 17, 1179–1196.

PubMed Abstract | Google Scholar

Crunelli, V., and Hughes, S. W. (2010). The slow (< 1 Hz) rhythm of non-REM sleep: a dialogue between three cardinal oscillators. Nat. Neurosci. 13, 9–17. doi: 10.1038/nn.2445

PubMed Abstract | CrossRef Full Text | Google Scholar

Cunningham, M. O., Whittington, M. A., Bibbig, A., Roopun, A., LeBeau, F. E., Vogt, A., et al. (2004). A role for fast rhythmic bursting neurons in cortical gamma oscillations in vitro. Proc. Natl. Acad. Sci. U.S.A. 101, 7152–7157. doi: 10.1073/pnas.0402060101

PubMed Abstract | CrossRef Full Text | Google Scholar

Dugladze, T., Schmitz, D., Whittington, M. A., Vida, I., and Gloveli, T. (2012). Segregation of axonal and somatic activity during fast network oscillations. Science 336, 1458–1461. doi: 10.1126/science.1222017

PubMed Abstract | CrossRef Full Text | Google Scholar

Fries, P. (2009). Neuronal gamma-band synchronization as a fundamental process in cortical computation. Annu. Rev. Neurosci. 32, 209–224. doi: 10.1146/annurev.neuro.051508.135603

PubMed Abstract | CrossRef Full Text | Google Scholar

Giraud, A.-L., and Poeppel, D. (2012). Cortical oscillations and speech processing: emerging computational principles and operations. Nat. Neurosci. 15, 511–517. doi: 10.1038/nn.3063

PubMed Abstract | CrossRef Full Text | Google Scholar

Goldental, A., Guberman, S., Vardi, R., and Kanter, I. (2014). A computational paradigm for dynamic logic-gates in neuronal activity. Front. Comput. Neurosci. 8:52. doi: 10.3389/fncom.2014.00052

PubMed Abstract | CrossRef Full Text | Google Scholar

Gray, C. M. (1994). Synchronous oscillations in neuronal systems: mechanisms and functions. J. Comput. Neurosci. 1, 11–38. doi: 10.1007/BF00962716

PubMed Abstract | CrossRef Full Text | Google Scholar

Grillner, S., Markram, H., De Schutter, E., Silberberg, G., and LeBeau, F. E. (2005). Microcircuits in action–from CPGs to neocortex. Trends Neurosci. 28, 525–533. doi: 10.1016/j.tins.2005.08.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Hasselmo, M. E. (2005). What is the function of hippocampal theta rhythm?—Linking behavioral data to phasic properties of field potential and unit recording data. Hippocampus 15, 936–949. doi: 10.1002/hipo.20116

PubMed Abstract | CrossRef Full Text | Google Scholar

Jirsa, V. K., and Haken, H. (1996). Field theory of electromagnetic brain activity. Phys. Rev. Lett. 77, 960. doi: 10.1103/PhysRevLett.77.960

PubMed Abstract | CrossRef Full Text | Google Scholar

Kahana, M. J. (2006). The cognitive correlates of human brain oscillations. J. Neurosci. 26, 1669–1672. doi: 10.1523/JNEUROSCI.3737-05c.2006

PubMed Abstract | CrossRef Full Text | Google Scholar

Kanter, I., Kopelowitz, E., Vardi, R., Zigzag, M., Kinzel, W., Abeles, M., et al. (2011). Nonlocal mechanism for cluster synchronization in neural circuits. EPL 93:66001. doi: 10.1209/0295-5075/93/66001

CrossRef Full Text | Google Scholar

Klimesch, W. (1996). Memory processes, brain oscillations and EEG synchronization. Int. J. Psychophysiol. 24, 61–100. doi: 10.1016/S0167-8760(96)00057-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Klimesch, W. (1999). EEG alpha and theta oscillations reflect cognitive and memory performance: a review and analysis. Brain Res. Rev. 29, 169–195. doi: 10.1016/S0165-0173(98)00056-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Marmari, H., Vardi, R., and Kanter, I. (2014). Chaotic and non-chaotic phases in experimental responses of a single neuron. EPL 106:46002. doi: 10.1209/0295-5075/106/46002

CrossRef Full Text | Google Scholar

Minlebaev, M., Colonnese, M., Tsintsadze, T., Sirota, A., and Khazipov, R. (2011). Early gamma oscillations synchronize developing thalamus and cortex. Science 334, 226–229. doi: 10.1126/science.1210574

PubMed Abstract | CrossRef Full Text | Google Scholar

Nir, Y., Mukamel, R., Dinstein, I., Privman, E., Harel, M., Fisch, L., et al. (2008). Interhemispheric correlations of slow spontaneous neuronal fluctuations revealed in human sensory cortex. Nat. Neurosci. 11, 1100–1108. doi: 10.1038/nn.2177

PubMed Abstract | CrossRef Full Text | Google Scholar

Roxin, A., Riecke, H., and Solla, S. A. (2004). Self-sustained activity in a small-world network of excitable neurons. Phys. Rev. Lett. 92:198101. doi: 10.1103/PhysRevLett.92.198101

PubMed Abstract | CrossRef Full Text | Google Scholar

Sanchez-Vives, M. V., and McCormick, D. A. (2000). Cellular and network mechanisms of rhythmic recurrent activity in neocortex. Nat. Neurosci. 3, 1027–1034. doi: 10.1038/79848

PubMed Abstract | CrossRef Full Text | Google Scholar

Silva, L. R., Amitai, Y., and Connors, B. W. (1991). Intrinsic oscillations of neocortex generated by layer 5 pyramidal neurons. Science 251, 432–435. doi: 10.1126/science.1824881

PubMed Abstract | CrossRef Full Text | Google Scholar

Thivierge, J.-P., Comas, R., and Longtin, A. (2014). Attractor dynamics in local neuronal networks. Front. Neural Circuits 8:22. doi: 10.3389/fncir.2014.00022

PubMed Abstract | CrossRef Full Text | Google Scholar

Vardi, R., Goldental, A., Guberman, S., Kalmanovich, A., Marmari, H., and Kanter, I. (2013a). Sudden synchrony leaps accompanied by frequency multiplications in neuronal activity. Front. Neural Circuits 7:176. doi: 10.3389/fncir.2013.00176

PubMed Abstract | CrossRef Full Text | Google Scholar

Vardi, R., Goldental, A., Marmari, H., Brama, H., Stern, E. A., Sardi, S., et al. (2015). Neuronal response impedance mechanism implementing cooperative networks with low firing rates and microsecond precision. Front. Neural Circuits 9:29. doi: 10.3389/fncir.2015.00029

CrossRef Full Text | Google Scholar

Vardi, R., Guberman, S., Goldental, A., and Kanter, I. (2013b). An experimental evidence-based computational paradigm for new logic-gates in neuronal activity. EPL 103:66001. doi: 10.1209/0295-5075/103/66001

CrossRef Full Text | Google Scholar

Vardi, R., Timor, R., Marom, S., Abeles, M., and Kanter, I. (2012a). Synchronization with mismatched synaptic delays: a unique role of elastic neuronal latency. EPL 100:48003. doi: 10.1209/0295-5075/100/48003

PubMed Abstract | CrossRef Full Text | Google Scholar

Vardi, R., Wallach, A., Kopelowitz, E., Abeles, M., Marom, S., and Kanter, I. (2012b). Synthetic reverberating activity patterns embedded in networks of cortical neurons. EPL 97:66002. doi: 10.1209/0295-5075/97/66002

CrossRef Full Text | Google Scholar

Wang, X.-J. (2010). Neurophysiological and computational principles of cortical rhythms in cognition. Physiol. Rev. 90, 1195–1268. doi: 10.1152/physrev.00035.2008

PubMed Abstract | CrossRef Full Text | Google Scholar

Wiest, M. C., and Nicolelis, M. A. (2003). Behavioral detection of tactile stimuli during 7–12 Hz cortical oscillations in awake rats. Nat. Neurosci. 6, 913–914. doi: 10.1038/nn1107

PubMed Abstract | CrossRef Full Text | Google Scholar

Wilson, H. R., and Cowan, J. D. (1972). Excitatory and inhibitory interactions in localized populations of model neurons. Biophys. J. 12, 1. doi: 10.1016/S0006-3495(72)86068-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, J.-Y., Guan, L., and Tsau, Y. (1999). Propagating activation during oscillations and evoked responses in neocortical slices. J. Neurosci. 19, 5005–5015.

PubMed Abstract | Google Scholar