Broadband macroscopic cortical oscillations emerge from intrinsic neuronal response failures

Broadband spontaneous macroscopic neural oscillations are rhythmic cortical firing which were extensively examined during the last century, however, their possible origination is still controversial. In this work we show how macroscopic oscillations emerge in solely excitatory random networks and without topological constraints. We experimentally and theoretically show that these oscillations stem from the counterintuitive underlying mechanism—the intrinsic stochastic neuronal response failures (NRFs). These NRFs, which are characterized by short-term memory, lead to cooperation among neurons, resulting in sub- or several- Hertz macroscopic oscillations which coexist with high frequency gamma oscillations. A quantitative interplay between the statistical network properties and the emerging oscillations is supported by simulations of large networks based on single-neuron in-vitro experiments and a Langevin equation describing the network dynamics. Results call for the examination of these oscillations in the presence of inhibition and external drives.

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 invitro 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.

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 substrateintegrated 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 2 . 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 2 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-5phosphonovaleric 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 60electrode 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 suprathreshold 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 firing frequency (i) = stimulation frequency · i+125 m=max(0,i−125) where Is_Spiked(m) = 1 if the neuron responded to stimulation no. m, otherwise Is_Spiked(m) = 0.

Simulations
Simulations (similar to Vardi et al., 2015) consist of a network of N leaky integrate and fire neurons 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-stimulationinterval, p(s|τ), should lead to <ISI> = τ c (Figures 1, 3). One can show that any p(s|τ) satisfying where p(τ) is the probability density of an inter-stimulationinterval 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 interstimulation-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 For this p(τ) some of the p(s|τ) solutions, which lead to <ISI> = τ c , are p(s|τ) = 0.5, p(s|τ) = τ/τ c and 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: 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.

Neuronal Plasticity: Neuronal Response Failures
We start with the quantification of the NRL (Vardi et al., 2013a(Vardi et al., , 2015Marmari et al., 2014), measured as the timelag 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 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  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.
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-timelags ( 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-stimulationintervals 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, τ.

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 interstimulation-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 subthreshold 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).

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., 2008Bollimunta et al., , 2011Nir 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 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.
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) −1 (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 selffeedback 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 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 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 realtime interface; AG performed the simulations and developed the theoretical framework with help of PS.; IK, RV, and AG wrote the manuscript.