Neuronal response impedance mechanism implementing cooperative networks with low firing rates and μs precision

Realizations of low firing rates in neural networks usually require globally balanced distributions among excitatory and inhibitory links, while feasibility of temporal coding is limited by neuronal millisecond precision. We show that cooperation, governing global network features, emerges through nodal properties, as opposed to link distributions. Using in vitro and in vivo experiments we demonstrate microsecond precision of neuronal response timings under low stimulation frequencies, whereas moderate frequencies result in a chaotic neuronal phase characterized by degraded precision. Above a critical stimulation frequency, which varies among neurons, response failures were found to emerge stochastically such that the neuron functions as a low pass filter, saturating the average inter-spike-interval. This intrinsic neuronal response impedance mechanism leads to cooperation on a network level, such that firing rates are suppressed toward the lowest neuronal critical frequency simultaneously with neuronal microsecond precision. Our findings open up opportunities of controlling global features of network dynamics through few nodes with extreme properties.


Introduction
The attempt to understand high cognitive functionalities and cooperative activities of neurons within a network results in many open questions. One question is which mechanism underlies the extremely low firing rates, few Hertz, of neurons comprising a network of threshold units, as a single neuron is capable of much higher firing frequencies (Amit and Brunel, 1997;Shafi et al., 2007;He et al., 2010;O'Connor et al., 2010). The second puzzle is the realization of precise neuronal response timings on a network level (Vanrullen et al., 2005;Butts et al., 2007;Panzeri et al., 2010), whereas their variations are typically in the order of several milliseconds (Lass and Abeles, 1975;Mainen and Sejnowski, 1995;Schoppa and Westbrook, 1999;Foust et al., 2010). This realization is a prerequisite for the feasibility of temporal codes, which might play a role in brain functionalities. In this work we demonstrate that network low firing rates and neuronal precise response timings both stem from a single neuronal property, the neuronal response impedance mechanism, which counter-intuitively leads to cooperation among individual neurons comprising a network.
The global features of network dynamics are governed by the distribution tail of the nodal properties, neuronal properties, rather than by specific distributions of the network links, the synapses.
The phenomenon of extremely low firing rates on the network level, about 1-3 Hz (Shafi et al., 2007;O'Connor et al., 2010), requires a balance between two opposing trends, spike birth and death. The rate of evoked spikes in an excitatory network is expected to constantly increase if each neuron excites several other neurons. This spike birth trend is moderated by a death trend, eliminating spikes by inhibitory synapses or weakening excitatory synapses to be sub-threshold (Turrigiano, 2008;Daqing et al., 2011). In case the birth and death trends are not precisely balanced, the network firing rate either diverges toward extremely high firing frequencies or practically vanishes (Turrigiano and Nelson, 2004;Chih et al., 2005). This balance can be theoretically achieved by several predefined synaptic designs depending on the topology of the network, such as a wide distribution of excitatory synaptic strengths balanced by a fraction of inhibitory synapses (Van Vreeswijk and Sompolinsky, 1998;Brunel, 2000;Song et al., 2000;Vogels and Abbott, 2009;Vogels et al., 2011;Teramae et al., 2012;Spiegel et al., 2014). A second approach relies on modular networks (Diesmann et al., 1999;Litvak et al., 2003;Kumar et al., 2010;Rad et al., 2012). The firing activity of each module is maintained by intra-excitatory synapses, whereas low firing rates are achieved by synaptic inhibition among the modules, resulting in their alternate firing.
The second phenomenon, the realization of consistent and reliable neuronal response timings, which is a prerequisite for a possible realization and usefulness of temporal coding, is limited by the least precise building block of the network, i.e., neurons and synapses. Synaptic conductance is reproducible with accuracy of dozens of microseconds (Csicsvari et al., 1998;Doyle and Andresen, 2001;Rodríguez-Moreno et al., 2011), which is meaningful for temporal coding only if the temporal precision of a neuron is comparable. Indeed, there is some experimental verification of specific neurons functioning with very fine temporal resolution (Carr, 1993;Agmon-Snir et al., 1998;Kayser et al., 2010); nevertheless, the temporal resolution of cortical neurons is unclear. Neuronal precision can be experimentally determined by the neuronal response latency (NRL), which reflects the internal dynamics of the neuron and is measured as the time-lag between a stimulation and its corresponding evoked spike (Wagenaar et al., 2004;De Col et al., 2008;Gal et al., 2010;Vardi et al., 2012). This quantity dynamically varies by several milliseconds and is expected to accumulate to dozens of milliseconds in a neuronal chain and even more in recurrent networks (Vardi et al., 2013a,b). Hence, a necessary prerequisite for the realization of consistent and precise temporal coding on a network level is in question. Nevertheless, the possible relevance and usefulness of temporal coding to brain functionalities are beyond the scope of our work.
In this work we study experimentally, in vitro and in vivo, the neuronal response impedance, i.e., the statistics of the neuronal response timings and response failures to different stimulation frequencies. Specifically we find that under low stimulation rates neuronal responses are precise up to several µs, while higher stimulation rates lead to an increased imprecision and response failures. Simulations of large networks as well as theoretical arguments supported by experimental evidences indicate that those neuronal response failures lead to low firing rates, even in excitatory networks. Hence, the dynamics on a network level lead to the coexistence of low firing rates and µs precision. It calls for the reexamination of the role of inhibition in the brain activity and its contribution to the stabilization of the network low firing rates. This work also implies that it is highly possible that neurons operate under µs precision, which is a prerequisite for the formation of reliable temporal coding.

Animals
All procedures were 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 (Vardi et al., 2012(Vardi et al., , 2013c. The cortical tissue was digested enzymatically with 0.05% trypsin solution in phosphatebuffered 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-4 weeks in vitro, prior to the experiments. 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 (i.e., ∼0.32 neurons in 100 µm 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 in vitro experiments, except for the experiments shown in Figure 9, 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 either 200 or 500 µm from each other (Multi-Channel Systems, Reutlingen, Germany) were used. The insulation layer (silicon nitride) was pre-treated with polyethyleneimine (0.01% in 0.1 M Borate buffer solution). A commercial setup (MEA2100-2x60-headstage, MEA2100interface 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 3channel 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,200] µ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 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 a duration of 200 µs using 30 repetitions at a rate of 5 Hz followed by 1200 repetitions at a rate of 10 Hz.

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.
Action potentials at experiments of the real-time adaptive algorithm for the stabilization of the neuronal response latency around a predefined latency were detected on-line by threshold crossing, using a detection window of typically 2-15 ms following the beginning of an electrical stimulation (Wagenaar et al., 2004).
In order to overcome the temporal precision of 20 µs determined by the maximal sampling rate of our recording device (50 kHz), in all in vitro experiments, except the above-mentioned adaptive algorithm, the following linear interpolation method for spike detection was used. For a given threshold crossing (V threshold ), we identify the two nearby sampling points: (t 1 , V 1 ) and (t 2 , V 2 ) where V 1 ≥V threshold , V 2 <V threshold and t 2 = t 1 + 20 µs. Using linear interpolation between these two sampling points, the threshold crossing time, t threshold , is estimated as (see also Figure 6) The neuronal response latency was then calculated as the duration from the beginning of a stimulation to t threshold .

Methods of Simulation
We simulated a network of 2000 excitatory leaky integrate and fire neurons (N = 2000). The voltage V' i (t) of neuron i (i ∈ [1,N]), is given by the equation: with V stable = −70 mV and a threshold of V threshold = −54 mV, τ = 20 ms is the membrane time constant, J ′ ji is the connection strength from neuron j to neuron i (see connectivity section) and D ji is the time delay from an evoked spike of neuron j to the stimulation of neuron i and is randomly chosen from a flat distribution U(6, 9.5) ms. Under the variable substitution V ′ = (V threshold − V stable )·V+V stable , J ′ ji = 16J ji mV and J ′ = 16J mV the equation of the voltage, V, becomes now For simplicity, we use this version of equation in the manuscript, since under this scaling V threshold = 1, V stable = 0, J > 1 is above threshold and J < 1 is bellow threshold. Nevertheless, results are the same for both equations. The initial voltage is V i (t = 0) = 0.5 ∀ i and the integration is done using the Euler method with 0.05 ms time step.
If V i crosses the threshold, 1, the neuron may fire (see response failure section). If the neuron fires the voltage is reset to −0.5 after a refractory period of 2 ms, in which the neuron is inactive, does not respond to new stimulations. In the case of a response failure, the voltage is set to 0.2 without a refractory period.

Connectivity
Initially all connections are set to zero, i.e., all J ji = 0. First, we go over all neurons and randomly select a post-synaptic neuron for each one. Each neuron can be selected as a postsynaptic neuron only once, and a neuron cannot be connected to itself, i.e., J ii = 0. As a result of this procedure each neuron has only one pre-synaptic neuron and only one post-synaptic neuron. Next, we select with a probability of 0.1/N, unless stated differently, additional above-threshold connections. The strength of all above-threshold connections is set to J ji = 2, i.e., J ′ ji =32 mV.

Response Failure
We define t(i,n) as the time neuron i crossed the threshold for the n th time. The response failure probability for the n th threshold crossing of neuron i is: where α is a measure of the neuronal forgetfulness and is equal to 1.4, unless stated differently. 1/τ C (i) is the critical frequency, f C (i), of neuron i, randomly chosen for each neuron, described in the figures. Note that negative response failure is taken as zero and in the limiting case of a periodic stimulation pattern, with the frequency f stim , the expected response failure is obtained, independent of n. The simulations on a network level are found to be independent of the initial conditions, which in the presented simulations were taken as P fail (i,0) = 0.

Initial Stimulations
To start the activity of the network the following external stimulations are chosen: First we randomly generate the external stimulation times for each neuron through a Poisson process with a rate of 50 Hz. Each external stimulation has a survival probability of exp (−T/200), where T stands for the time of the external stimulation, measured in milliseconds. Only stimulations with T < 1000 ms are taken into account, as after 1000 ms the network has a self-sustaining activity, without spontaneous firing. Results are found to be insensitive to different initial conditions, e.g., only a single stimulation to a single neuron at time 0.

Frequency Histograms
The firing frequency of a neuron was determined as the number of times the neuron fired between 4 and 59 s divided by 55 s. The bin size for all histograms is 0.5 Hz ( Figures 11B-E, 12).

Additional Response Failure Probability
In the simulations shown at Figure 11D the response failure probability is C+(1-C) ·P fail , where P fail is given by Equation (1) Frontiers in Neural Circuits | www.frontiersin.org and C is a constant. For stimulation frequencies below f C , the response failure probability is C as P fail is zero. In Figure 11D, C = 0.07.

Spontaneous Activity
In the simulations shown at Figure 11E the probability for a spontaneous stimulation per integration time step is 5·10 −5 , resulting on the average in 1 spontaneous stimulation per second for each one of the neurons.

Additional Sub-Threshold Connections
In the simulations shown at Figure 12C there are additional subthreshold connections, where every J ji = 0 was changed to J ji = J sub with a probability P sub (J ii = 0 ∀ i).

Surgery
The experiments were performed on Sprague-Dawley or Wistar rats weighing 100-200 g, initially anesthetized with urethane (1.25 g per kg body weight, intraperitoneal), and given hourly supplemental injections of ketaminexylazine (30 and 7 mg per kg, respectively, intramuscular). The animals' temperatures were maintained at 37 • C, while placed in a stereotaxic instrument. All animals continued to breathe without artificial respiration, and were suspended by ear-bars and a clamp at the base of the tail to minimize movements of the brain caused by breathing. The cisterna magna was opened to relieve intracranial pressure. Small (1 mm) holes were drilled in the skull to allow insertion of electrocorticogram electrodes, in addition to a 8 × 8 mm craniotomy for insertion of the recording and stimulating electrodes over the prefrontal cortex, where the dura was removed.

Electrophysiology
Two stimulating electrodes (tapered tungsten wire) were inserted about 0.6 mm deep into the cortex, to approximately cortical layer IV. Their tips had fixed distances of about 0.2 mm. The intra-cellular recording electrode was inserted within a radius of 1.5 mm into a depth of 0.2-1.2 mm, as a glass micropipette filled with 2-4% biocytin (Sigma) dissolved in 1 M potassium acetate. Electrode resistances ranged from 30 to 100 M . After insertion of both types of electrodes, the exposed cortex was covered with a low-melting-point paraffin wax to reduce brain pulsations. Stimulus current amplitudes used were between 50 and 700 µA, and the duration of current flow was 100-200 µs. Recordings were made using an active bridge amplifier and then filtered and digitized at a rate of 10 kHz. The experiments were performed on neurons with membrane potentials lower than −60 mV and action potentials higher than 0 mV. at stimulations (2n + 1)·100 τ =197 ms and at (2n + 2)·100 τ =37 ms, where n = 0,1,.... Results indicate an asymmetry between facilitation, with an amplitude of ∼0.4 ms in case τ decreased, and depression, with an amplitude of ∼0.1 ms in case τ increased (see zoom-in, gray). (D) Same as (C) but with moderate changes in τ, at stimulations (2n + 1)·100 τ = 150 ms and at (2n + 2)·100 τ = 84 ms. Results indicate an asymmetry between facilitation, with an amplitude of ∼0.2 ms in case τ decreased, and depression, with an amplitude of ∼0.1 ms in case τ increased (see zoom-in, gray).

Stabilization of the Neuronal Response Latency
When a neuron is repeatedly stimulated its response latency, NRL, stretches gradually. This effect was demonstrated for a 30 Hz stimulation frequency ( Figure 1A) using cultured cortical neurons functionally separated from their network by synaptic blockers and stimulated such that most stimulations resulted in an evoked spike (Materials and Methods). The accumulated stretching over few hundreds of repeated stimulations is typically several milliseconds, comparable with the initial NRL ( Figure 1A). This stretching terminates at the intermittent phase, where the average NRL remains constant and is accompanied by both large fluctuations that can exceed a millisecond and a non-negligible fraction of neuronal response failures ( Figure 1A). The NRL increase is a fully reversible phenomenon , which considerably decays after few seconds without stimulations.  The electrode-amplified noise, σ, is presented by the horizontal bars around the blue circle. This uncertainty in the electrode voltage induces an uncertainty at the interpolated threshold crossing time, t, following the relation t σ = 20 V . The time-lag between two samples is 20 µs (50 kHz sampling rate), and the average voltage difference between two samples in this spike phase is ∼15 µV. For the above-mentioned typical values, σ = 2.14 µV and V = 15 µV, one finds t = 2.85 µs.
Frontiers in Neural Circuits | www.frontiersin.org The stretching of the NRL seems to prevent consistent response timings, required for the realization of temporal codes. Even at the intermittent phase, where the average NRL is stable, the emergence of large fluctuations and neuronal response failures opposes the realization of temporal codes. Hence, we examine the feasibility of NRL stabilization using three different stimulation scenarios.
The proof of concept for NRL stabilization is examined first under an adaptive external feedback scenario (Steingrube et al., 2010) (Figure 1B). It relies on the fact that the average stretching of the NRL per stimulation increases with the decrease of the time-lag between stimulations, and vice versa (Vardi et al., 2012). Accordingly, a real-time adaptive algorithm for the stabilization of the NRL around a predefined latency, L ST , was experimentally tested, indicating a supreme stabilization measured by the standard deviation σ ∼16 µs (Figure 1C). In order to maintain stabilization around L ST the time-lag between stimulations has to be adjusted continuously, however, after ∼1000 stimulations the average τ stabilizes around 200 ms.
This proof of concept for NRL stabilization, without response failures, calls for a more natural realization of this phenomenon without adaptive external control. Using a fixed time-lag of ∼200 ms between stimulations, equal to the average timelag obtained under the external feedback algorithm (Figure 1C), stabilization around the same L ST ∼5 ms is obtained, with a comparable standard deviation of σ ∼14 µs (Figure 1C), however the transient to stabilization is longer (Figures 2A,B). Results suggest that the stabilized NRL, L ST , is primarily a function of the average time-lag between stimulations. This was confirmed through a more realistic scenario where the same neuron was stimulated following random time-lags characterized by the same average of ∼200 ms (Figure 1C), resulting in a similar L ST with slightly larger fluctuations, σ. In this case of stimulation with random time-lags, the L ST is slightly reduced as a result of asymmetric fluctuations in the neuronal response to sudden decrease or increase of the time-lags between consecutive stimulations, forming momentarily depression or facilitations in the neuronal response latency (Figures 2C,D).
We examine the effect of an abrupt transition to a much higher stimulation frequency following different periods of stabilization around L ST . This scenario is exemplified by the stabilization of the neuronal response latency around L ST ∼9.7 ms using τ = 110 ms, followed by an abrupt change to much shorter timelags between stimulations, τ = 25 ms ( Figure 3A). A comparison between the segments of the neuronal response latency profiles following the change in the stimulation frequency shows that they are fairly identical (Figure 3B), indicating a lack of dependency on the history of simulations. In contrast, when the same neuron is stimulated solely with τ = 25 ms its neuronal response latency profile above L ST is clearly different (Figures 3A,B). Before the Frontiers in Neural Circuits | www.frontiersin.org latency stabilizes, it is solely a function of a global quantity, the averaged time-lag between all previous stimulations constitute the current neuronal response latency stretching . This behavior indicates neuronal long-term memory , in contrast to the behavior after stabilization.
The silencing of long-term memory during the stabilized period of the neuronal response latency does not disable the ultra-fast neuronal plasticity. Specifically, a substantial decrease in τ results in a momentary decrease in the neuronal response latency, indicating ultra-fast plasticity in the form of facilitation (Figures 3A,C). For example, when τ is shortened from 110 to 25 ms, neuronal facilitation with an amplitude of ∼0.2 ms is evident, and the subsequent latency profiles are independent of the duration of stabilization at L ST (Figure 3C). A similar facilitation, with an amplitude of ∼0.2 ms, also occurs when τ is momentarily shortened from 110 to 25 ms once every 100 stimulations (Figure 3D). The continuation of the response latency profile is recovered immediately following the next stimulation, either during the latency stretching phase or at the stabilization around L ST , indicating insensitivity of global neuronal response latency profiles to momentary leaps in τ. The co-existence of these two features, silenced long-term neuronal memory and neuronal ultra-fast plasticity (Figure 3A), suggests a possible realization of reliable signaling in temporal codes , especially when the neuronal response latency is stable.

In Vitro Experiments
The NRL profiles differ among neurons in their total increase, average increase per stimulation and detailed profile forms. In case the NRL can be stabilized at any latency, are there common neuronal trends, independent of the detailed NRL profile? Specifically, can one find universal features characterizing the fluctuations around L ST , σ(L ST ), and the average time-lag between stimulations required for stabilization around L ST , τ(L ST )? Stabilization close to the initial NRL requires τ in the order of seconds which rapidly decreases to typically ∼[50,150] ms while approaching the NRL at the intermittent phase (Figures 4, 5A,D). These two regions are also reflected in the behavior of σ(L ST ), where below a certain NRL, approximated by L T , the standard deviation, σ, is almost constant and diverges as a power-law as the intermittent NRL is approached (Figures 4, 5C,F). The maximal σ varies much among neurons (Figures 5B,E), and similarly large deviations in the exponent characterizing the power-law are observed (Figures 5C,F). Nevertheless, the power-law was found to be stimulations (x-axis) for the neuron in (C) stimulated at 12 Hz (orange dots), where 1/0 stand for evoked spike/response failure. These probabilities are compared to the corresponding theoretical probabilities (black triangles) assuming uncorrelated failures with the measured P fail ∼0.4. All measurements were taken at the intermittent phase. (E) The probabilities for the occurrence of segments of m consecutive response failures bounded by evoked spikes, P 0 (m), measured at the intermittent phase for the neuron in (C) stimulated at 12 Hz (orange dots). The theoretical values, ln(1-P fail )+m·ln(P fail ) with P fail = 0.4 (black line), similar to (C). (F) Similar to (E), the probabilities for the occurrence of segments of m consecutive evoked spikes bounded by response failures, P 1 (m), measured at the intermittent phase for the neuron in (C) stimulated at 9 Hz (orange dots). The theoretical values, ln(P fail )+m·ln(1-P fail ) with P fail = 0.15 (black line), similar to (C). All experiments shown in this figure were done on a cultured neuron.
Frontiers in Neural Circuits | www.frontiersin.org a universal characteristic of the transition to the intermittent phase.
The NRL at the transition between the regions of semiconstant and rapidly increasing σ, L T , is estimated where the NRL profile changes its concave form to convex, d 2 L/dStimulation 2 = 0, under relatively high stimulation frequencies (Figures 5A,D). The concave profile, L<L T , was recently found to identify a reproducible non-chaotic neuronal phase . In this phase deviations between NRL profiles, obtained in different stimulation trials consisting of identical stimulation pattern of a single neuron, do not increase with the number of stimulations . In an opposite manner, for a convex NRL profile, L>L T , a chaotic neuronal phase emerges, indicating an exponential divergence among NRL profiles obtained in different trials. The emergence of nonchaotic and chaotic neuronal phases preceding the intermittent phase  is also reflected in the NRL stabilization. Stabilization at L ST <L T is realized with semiconstant small σ, whereas at L ST >L T a rapid increase in σ is observed.
The minimal standard deviation of few microseconds, i.e., σ∼5 µs (Figures 5B,E), is measured when L ST is close to the initial NRL, and requires sub-Hertz stimulation rate. These minimal deviations saturate our experimental lower bound of σ, which stems from the unavoidable amplified noise measured by the electrode (Figure 6), thus the examination of a better neuronal stabilization is beyond our experimental limitations. To overcome these limitations, we turn to the following theoretical argument. The number of ions in an evoked spike can be estimated from the density of the neuronal membrane capacitance which is around 1 µF/cm 2 (Cole, 1968;Gentet et al., 2000). Hence, for the diameter of a typical neuronal soma, ∼20 µm (its surface area ∼10 3 µm 2 ) (Lübke et al., 1996), one finds the neuronal membrane capacitance as C∼10 −11 F. Using the voltage difference in the membrane during an evoked spike (Cole, 1968), V∼0.1 V, the scaling of the number of ions in an evoked spike is estimated as Q = C V∼10 7 e. These ∼10 7 ions are evoked in about 1 ms, the duration of a spike, implying that an individual ion is evoked on the average every 1 ms/10 7 =10 −10 s. Assuming a simple stochastic process for the sequential emission of ions, the expected time deviation for an evoked spike is t∼10 −10 √ 10 7 ∼0.3 µs, which is only slightly below the experimentally demonstrated σ of extremely few microseconds.

In Vivo Experiments
The NRL stabilization, with fluctuations in the order of 100 µs, was confirmed also in in vivo experiments (Figure 7). Results The experimentally measured probability for response failures, P fail (orange dots) and the corresponding τ (crimson dots) at different frequencies in the range of [6,30] Hz. For each frequency the averaged ISI= τ/(1-P fail ) (crimson circles) is very close to τ C = 160 ms (f C ∼6.2 Hz), independent of the stimulation frequency. support the in vitro experiments such that σ can reach ∼100 µs without response failures despite 1-8 ms NRL stretching. However, in enhanced stimulation frequencies the intermittent phase emerged, characterized by increased fluctuations and the FIGURE 10 | Neuronal impedance mechanism on a chain level leading to low firing rates. (A) Schematic (top) of a cultured neuron with f C = 6.3 Hz and its stimulation/firing frequencies, f1/f2, respectively. Stimulation frequency (red) and firing frequency (green) smoothed using 100 Stimulation sliding window of a neuron stimulated with alternating time-lags, 20 and 300 ms, equal on the average to ∼τ C , each time-lag repeats 5 to 10 times before switching. The averaged f1 = 6.2 Hz (red dashed-line) and the theoretically predicted f2 = 4.3 ms (green dashed-line) are shown for comparison. The resulting firing frequency of the cultured neuron is ∼4 Hz (green full line), substantially lower than its f C , and is close to the predicted f2 following ISI = 0.5(160 + 300) = 230 ms. (B) A cultured neuron with f C ∼5.5 Hz was stimulated periodically, 40 times at 12 Hz and 40 times at 7 Hz (inset). The probability of a response failure, P fail (i), i = 0, …,79 (orange dots), measured over 200 recurrences, fitted with optimized α = 1.1 following equation (3) (black line) and the expected standard deviation P fail (i) 1 − P fail (i) /200 (gray area). The experimentally observed probabilities are similar to the predicted values under a fixed stimulation frequency, P fail = (τ C − τ)/τ C = 1-f C /f, resulting in 1-(5.5/12)∼0.54 for 12 Hz and 1-(5.5/7)∼0.21 for 7 Hz.  , of the firing rates were found to be independent of the initial external stimulation patterns given to the network, indicating neuronal cooperation that reduces firing frequencies toward the lowest critical frequencies, f C . All the results shown in this figure were produced in simulations.
emergence of response failures. A precise characterization of the intermittent phase requires much longer stimulation periods, which are limited in our experimental setup, and also taking into account the neuronal spontaneous activity, which in the presented neuron was ∼2 Hz. The ∼1.5 mm distance between the stimulating and recording electrodes (Abeles, 1991;Kincaid et al., 1998;Zheng and Wilson, 2002) and the large (∼17 ms) NRL stretching at the intermittent phase ( Figure 7D) support a synaptic mechanism with high probability (through orthodromic stimulation) and probably even relaying via several synapses (Abeles, 1991;Kincaid et al., 1998;Zheng and Wilson, 2002). Note that the lack of a clear L T might be attributed to the accumulation of the NRL along a neuronal chain. These finding strongly support the temporal robustness of signal transmission in the brain (Csicsvari et al., 1998;Bonifazi et al., 2005;Boudkkazi et al., 2011), even in cases where both synaptic and neuronal transmission are involved.

Universal Properties at the Intermittent Phase
Above a critical stimulation frequency, f C , the neuron enters the intermittent phase, characterized by a maximal average NRL, L C , independent of the stimulation frequency (Gal et al., 2010), accompanied by the appearance of response failures (Figures 1A,  8A) and large σ (Figures 5B,E). Specifically, there is a critical time-lag τ C (corresponding to f C ) where L C is first achieved accompanied by large σ but with a vanishing fraction of response failures ( Figure 8B). The transient times to reach stable NRLs are comparable for f≤f C (Figures 1C, 4, 8A,B). Nevertheless, as f C is approached the transient time to a stationary σ is much enhanced (Figure 8B), together with a power-law divergence of σ (Figures 5C,F), as expected in such a second order phase transition (Stanley, 1987). At stimulation frequencies above f C , neuronal response failures stochastically emerge such that their fraction is well approximated by P fail = (τ C − τ)/τ C (Figure 8C). Consequently, a global tenable quantity, the average inter-spike-interval (ISI), is preserved and is equal to τ/(1-P fail )∼τ C (Figure 8C), hence the neuron functions similar to a low pass filter. The stochastic occurrence of response failures was quantitatively examined using the following two tests. The first consists of calculating the probability for the occurrence of all combinations of responses for a given set of consecutive stimulations. These probabilities are then compared to the expected values under the assumption of random uncorrelated response failures with a given P fail (Figure 8D). The second test consists of calculating the probabilities, P 0 (m), for the occurrence of segments of m consecutive response failures bounded by evoked spikes. These probabilities were found to be in a good agreement with the theoretically predicted ones based on a Poissonian process with a rate -ln(P fail ) ( Figure 8E). Stimulation frequencies close to f C result in P fail <<0.5, and enable the examination of stochastic response failures for much longer segments consisting of m evoked spikes bounded by response failures, P 1 (m) (Figure 8F). Both tests indicate that response failures emerge stochastically and independently at the intermittent phase.
Similar results for the neuronal critical frequency, f C , as well as stabilization of the NRL were also found for cultured cortical neurons that were not functionally separated from their network by the addition of synaptic blockers (Figure 9), strengthening the biological relevance of our findings.

Fast Neuronal Adaptation to Frequency Modulation
The saturated neuronal firing frequency, f C , functions as an impedance mechanism limiting the average neuronal firing rate. Typically, f C is in the range of 6-15 Hz, but can be extended for some neurons as high as 27 Hz and as low as 3 Hz.
A plausible biological scenario that suppresses the firing frequency of a single neuron below f C is aperiodic time-lags between stimulations, as verified experimentally ( Figure 10A). For illustration, assume a slow mode of alternation between stimulation frequencies of 2f C (0.5τ C time-lag between stimulations) and 2f C /3 (1.5τ C ), such that the average time-lag between stimulations is τ C . For the high and low frequency modes, the expected probability for response failures is 0.5 and 0, respectively. Consequently, the average ISI is 0.5(1.5 τ C + τ C ) = 1.25τ C , corresponding to a lower firing rate of 0.8f C .
The fairly good agreement between the lowered firing rate, below f C , and the predicted one ( Figure 10A) strongly indicates fast neuronal plasticity (adaption) where the probability for a neuronal response failure is intrinsically adjusted following the temporary stimulation frequency. To quantify the time scale of this type of neuronal plasticity, a neuron characterized by f C ∼5.5 Hz was repeatedly stimulated with a recurrence of 80 stimulations, 40 at 12 Hz and 40 at 7 Hz ( Figure 10B, inset). The probability for response failure for each of the 80 stimulations, measured over 200 recurrences, consists of two semi-stationary values; ∼0.6 is attributed to stimulations given at 12 Hz and ∼0.3 to stimulations given at 7 Hz ( Figure 10B). The transient time between these values was found to vary between five and several dozen stimulations among neurons.
The probability profile of response failures, P fail , (Figure 10B) was fitted to the function where the stimulation number i is an integer in the range [0,79], τ (m) is the time-lag between stimulations m and m + 1, a negative m reads as τ (m+80) and the term (τ C − τ (m) )/τ C represents the probability for a response failure following τ (m) . The last term represents the weighted exponential decay function with the optimized fitted fading coefficient α ( Figure 10B) and A is the normalization coefficient setting P fail (i)=1 if all τ (m) are zero (see Materials and Methods). This quantitative modeling of intrinsic short-term plasticity enables the examination of its abundant cooperative effects within neuronal chains and networks.

Neuronal Impedance Mechanism on a Network Level
We experimentally examined a chain of two neurons, characterized by the critical frequencies f C = 7.5 and 15.5 Hz for the first and second neuron, respectively ( Figure 10C). The first neuron was stimulated with time-lags between stimulations taken randomly from U(20, 110) ms (Figure 10C), resulting in a ∼7.3 Hz average firing rate, close to its f C , and consequently its σ is large, exceeding 350 µs (Figure 10C). The second neuron is stimulated by the first neuron at ∼7.3 Hz on the average, far below its f C = 15.5 Hz. Hence, it relays the stimulations in the form of evoked spikes without response failures (P fail ∼0) and with supreme precision, σ ∼40 µs ( Figure 10C).
The effect of the neuronal response impedance mechanism on the network level was examined using large scale simulations of excitatory networks composed of N = 2000 leaky integrate and fire neurons (see Materials and Methods), whose prototypical topology was constructed using the following two steps. Each neuron is first randomly selected to have exactly one post-synaptic and one pre-synaptic connection, and the remaining connections are then selected with a survival probability of 0.1/N ( Figure 11A). All connections are abovethreshold, delays are taken randomly from U(6, 9.5) ms and neurons are selected randomly to have either f C ∼6.66 Hz (τ C = 150 ms) or f C ∼14.28 Hz (τ C = 70 ms) ( Figure 11B). Response failures for each neuron were implemented using the impedance mechanism, Equation (1) with α = 1.4 (see Materials and Methods), however results were found to be insensitive to the FIGURE 12 | Robustness of low firing rates on a network level.
(A-C) The normalized histogram of critical frequency, f C (light blue), and firing frequency (dark blue). (A) The low firing activity is robust to variations in α. An identical network topology as in Figure 11C but with different α: 1.4 (upper panel), 1.1 (middle panel) and 0.5 (lower panel), with mean firing rate of 6.26, 6.43, and 6.96 Hz, respectively. (B) The low firing activity is robust to higher connectivity. Networks obeying the same statistical features as in Figure 11C but the probability for an additional above-threshold connection, P, (see Materials and Methods) is 0.1/N (upper panel), 0.3/N (middle panel), and 0.5/N (lower panel), with mean firing rate of 6.26, 7.43, and 8.23 Hz, respectively. Results indicate that higher connectivity slightly increases the firing rates. (C) The low firing activity is robust to additional sparse/dense sub-threshold connectivity. Networks obeying the same statistical features as in Figure 11C where each neuron has, on the average, additional N·P sub =3 (left panel), 0.1·N (right panel) sub-threshold post-synaptic connections with the strength of J sub = 0.4 (left panel), 12/N (right panel). Note that the average total strength of the sub-threshold connections, per neuron, is preserved for both sparse and dense cases as N·P sub ·J sub = 1.2. All the results shown in this figure were produced in simulations.
FIGURE 13 | The effect of inhibition on the network spontaneous activity. The average spontaneous spike detection rate recorded from all 60 electrodes of a plated MEA over 10 min where no external stimulations were given (blue). Similarly, the spontaneous activity of the same MEA with addition Bicuculline (5 mM), which blocks inhibition (green). The average spike detection rates recorded from all electrodes (excluded the grounded 15th electrode) are presented by the dashed lines for each scenario.
precise α (Figure 12A). The distribution of firing rates of the 2000 neurons was estimated after several seconds using a time window of about 50 s (Figures 11B-E, 12), indicating ∼5.4 Hz averaged firing rate ( Figure 11B) which is even below the minimal f C = 6.66 Hz. This result is a consequence of the two abovementioned experimentally verified effects. The first, a single neuron along the chain, characterized by f C = 6.66 Hz, enforces all its consecutive neurons to fire no higher than this frequency. Since in a recurrent network most neurons have an ancestor neuron with such low f C = 6.66 Hz, they are expected to lower their firing rates toward this frequency. The second effect is aperiodic time-lags between stimulations, leading to firing rates even further below f C = 6.66 Hz ( Figure 11B).
The cooperative effect of the neuronal response impedance mechanism drives the average firing rate of the entire network even below the lower tail of the distribution of f C . In this state of low firing rates, neurons are typically in the non-chaotic phase, where deviations of only several microseconds around the average NRL are expected ( Figure 5). Hence, low firing rates and microsecond neuronal precision are simultaneously achieved on a network level even in the presence of neurons which can potentially fire at very high rates (McCormick et al., 1985;Tateno et al., 2004). This tendency was found to be robust to a more realistic scenario where f C was taken randomly from U(6.66, 14.28) Hz ( Figure 11C). Since our experiments indicate that f C can be as low as 3 Hz, compared to 6.66 Hz used in the above simulations, in large scale neural networks even lower firing rates are expected, as experimentally observed in cortical activity (Shafi et al., 2007;O'Connor et al., 2010).
The lack of response failures for stimulation frequencies below f C is too simplistic of an assumption for cortical dynamics, as failures may be generated, for instance, by synaptic noise and background inhibition. A theoretical argument and simulations ( Figure 11D) indicate that indeed the proposed cooperative effect is robust to an additional response failure probability, p < 0.075, for all frequencies, including those below f C . This critical probability, 0.075, is a result of the average chain length, which for our network topology is ∼9, as the probability for a neuron to have two post-synaptic connections is 0.1. Consequently, the probability for an evoked spike from the last neuron in the chain, given a stimulation to the first one is (1-p) 9 . Since the chain terminates in a branch to two consecutive chains, the preference of spike birth over spike death requires 2(1-p) 9 >1, resulting in p < 0.075. Similarly, low firing rates on a network level were found in simulations to be robust to the scenario of spontaneous stimulations at an average rate of 1 Hz per neuron (Figure 11E), as well as for the same network architecture with additional above-and sub-threshold connections, for both sparse and dense scenarios (Figures 12B,C).

Discussion
The dynamical properties of networks are typically assumed to reflect the statistical properties of their links. Following this framework, the low firing rate for a given neural network topology was achieved in simulations using specific distributions of excitatory and inhibitory synapses, which balance spike birth and death trends. Since an inhibitory synapse, a directed link, probabilistically blocks an evoked spike of its driven node only in a given time window (Vardi et al., 2013b), the low firing rates are expected to be sensitive to small changes in network topology, synaptic delays and to the emergence of spontaneous activity (Daqing et al., 2011), unless a mean-field limit is assumed. In this work we experimentally present the neuronal response impedance mechanism on a single neuron level, which results in stochastic neuronal response failures at high stimulation frequencies and in precise response timings at low stimulation frequencies. On a network level, this mechanism leads to robust low firing rates, where each node, neuron, independently generates response failures above a critical stimulation frequency and functions similar to a low pass filter. Consequently, cooperation among individual neurons, enforced by the network dynamics, results in low firing rates which are governed by the low critical frequencies of the extreme nodes. As a byproduct of these low firing rates, the nodal response timings are stabilized with microsecond neuronal precision. The emergence of the two cooperative features on a network level supports the possibility that all building blocks of neural networks, neurons and synapses, jointly operate under the same extreme precision.
An indirect experimental support that low firing rate on the network level is achieved using solely neuronal response failures and without inhibition is also presented (Figure 13). We compare the spontaneous activity of the same MEA plated with cultured cortical neurons without additional synaptic blocker and with additional Bicuculline which blocks inhibition. In both cases, the average spontaneous activity of each electrode is measure over 10 min. Results, exemplified in Figure 13, clearly indicate that in both cases low firing rates are maintained, although the suppression of the inhibition slightly enhanced the spike detection rates of most of the electrodes. It is clear that typically each electrode records the spiking activity of more than one neuron, hence the firing activity recorded by the MEA does not directly count the activity of the entire neural network. Nevertheless, the only slight increase in average firing activity measure by the MEA, supports our findings that the phenomenon of low firing rates is mainly attributed to the neuronal response failures.
The question arises as to what functionalities demand synaptic inhibition, inhibitory synapses. A possible hypothesis is that stationary network activity of low firing rates acts as a baseline cortical state. Over this state of activity, meaningful neuronal functionalities are embedded by the conditional temporal formation of neuronal firing, resulting in effective spatial summation, opening new routes of information flow through the network. Alternatively, temporary higher frequency stimulations to a subset of neurons result in abrupt changes of their NRLs (Figure 3). The accumulation of these changes along neuronal pathways (Vardi et al., 2013a(Vardi et al., ,b, 2014Goldental et al., 2014) dynamically changes the topography of the network (Vardi et al., 2013b;Goldental et al., 2014). Both abovementioned scenarios only temporarily influence the baseline cortical state, as the NRL stretching is a fully reversible phenomenon (De Col et al., 2008;Vardi et al., 2013c).
This work does not contradict the known mechanisms, e.g., inhibition, which lead to low firing rates. Nevertheless, the proposed mechanism is significantly robust to changes in the network structure, e.g., connectivity and synaptic strengths. Specifically, we have shown that neural network will exhibit low firing rates even without inhibition, as the proposed mechanism is based on the unreliable responses of the neurons. Our findings call for the reexamination of the role of inhibition as the main suppressor of firing rates in neural networks. Specifically, what is the synergism between inhibition and the intrinsic neuronal impedance mechanism and in what dynamical circumstances one is more dominant than the other? Furthermore, we have shown that although neurons generally do not respond in a temporal precision and imprecisions accumulate along the network, the low firing rates lead to a supreme stabilization of the neuronal responses. This stabilization is characterized by imprecisions of only several µs, opening the doors for the feasibility of temporal coding in various functionalities in the brain.

Author Contributions
RV and HM prepared the cultural tissues, the experimental materials and performed the in vitro experiments. AG performed the simulations and participated in the development of the theoretical arguments. RV, HM, and AG analyzed the data. SS participated in some of the in vitro experiments. PS designed the interface for the real-time experiments. HB, ES, and IK designed, performed and analyzed the in vivo experiments. RV and AG prepared the manuscript. IK initiated the study and supervised all aspects of the work. All authors discussed the results and commented on the manuscript.