Brief wide-field photostimuli evoke and modulate oscillatory reverberating activity in cortical networks

Cell assemblies manipulation by optogenetics is pivotal to advance neuroscience and neuroengineering. In in vivo applications, photostimulation often broadly addresses a population of cells simultaneously, leading to feed-forward and to reverberating responses in recurrent microcircuits. The former arise from direct activation of targets downstream, and are straightforward to interpret. The latter are consequence of feedback connectivity and may reflect a variety of time-scales and complex dynamical properties. We investigated wide-field photostimulation in cortical networks in vitro, employing substrate-integrated microelectrode arrays and long-term cultured neuronal networks. We characterized the effect of brief light pulses, while restricting the expression of channelrhodopsin to principal neurons. We evoked robust reverberating responses, oscillating in the physiological gamma frequency range, and found that such a frequency could be reliably manipulated varying the light pulse duration, not its intensity. By pharmacology, mathematical modelling, and intracellular recordings, we conclude that gamma oscillations likely emerge as in vivo from the excitatory-inhibitory interplay and that, unexpectedly, the light stimuli transiently facilitate excitatory synaptic transmission. Of relevance for in vitro models of (dys)functional cortical microcircuitry and in vivo manipulations of cell assemblies, we give for the first time evidence of network-level consequences of the alteration of synaptic physiology by optogenetics.

Since first expressing channelrhodopsin in neurons 1 , optogenetics 2 raised immense interest as one of the most promising techniques in neuroscience [3][4][5][6][7][8] . Optogenetics enables the manipulation of the activity of genetically-identified neurons with unprecedented temporal resolution [9][10][11][12][13][14] . Therefore, both in fundamental research and neuroprosthetics, the idea of regulating neuronal firing by light (e.g. in a closed loop) quickly emerged as a pivotal causative approach to advance our correlative understanding of brain (dys)functions. Nevertheless, the limits and most useful features of the stimulation (e.g. intensity, time course, frequency, temporal pattern, etc.), to be exploited in a control system, have not been exhaustively explored, but see 15,16 .
Despite the low spatial resolution in most in vivo applications, but see 17,18 , as light broadly addresses the feed-forward projections expressing an opsin, the target neuronal populations can be efficiently controlled 5,19,20 . In this case, optogenetic control of downstream neuronal firing is straightforward and irreplaceably elegant and simple. However, when recurrent collaterals expressing an opsin dominate, e.g. as in the cortex, it is expected that light stimuli elicit additional reverberating responses, due to local recurrent pathways. These responses arise as second-order effects from the recruitment or suppression of neuronal activity and may reflect unanticipated time-scales, and complex dynamical microcircuit properties. These add complexity to the interpretation of optogenetic stimulation and to its use for closed loop control, particularly in vivo where they might be hard to isolate and study in depth. A complementary approach for dissecting reverberating activity may come from using of reduced in vitro experimental preparations [21][22][23] , offering substantial advantages over in silico simulations. In fact,

Results
We cultured large-scale networks of rat primary cortical neurons, over 57 arrays of microelectrodes (MEAs) integrated in transparent glass substrates. Over four weeks in vitro, neurons developed functional synaptic connections and reached a mature pattern of spontaneous activity 34 . By detecting action potentials recorded at each microelectrode, we investigated the collective neuronal responses upon brief wide-field light stimulation (Fig. 1A). In 39 MEAs, we transduced cells by AAV to express a variant of channelrhodopsin (ChR2 LC-TC) fused to a red-fluorescent protein (mCherry), targeting selectively CaMKIIα -positive cells (i.e. putative glutamatergic neurons) 40 (Fig. 1D). This mutant opsin responds to blue light with stronger currents, reduced proton permeation, and enhanced calcium (Ca 2+ ) selectivity, compared to its wild type 41 . Immunocytochemistry confirmed our desired high level of opsin expression, transduction efficiency (i.e. ~80% of all neurons), and low apoptosis (not shown). Routine MEA recordings showed no differences with control cultures in terms of network development and cell survival, inferred by analyzing spontaneous bursting activity (Fig. 1B,C) 22,42 and by quantifying the number of active microelectrodes (i.e. those detecting > 1 spike in 50 s; Fig. 1C).

Light-evoked network responses.
We employed a computer-controlled blue LED to deliver brief wide-field pulses of light over the inner bottom area of each MEA. Simultaneously we recorded the neuronal spike responses, detected extracellularly at each of the 59 microelectrodes of the MEA (Fig. 1A, inset). The duration of light pulses was systematically varied in the set {0.1, 0.5, 1, 2, 5, 10, 20} ms, while the light intensity was kept constant to its full power value (i.e. 2 mW/mm 2 at the sample). This allowed us to vary reliably the amplitude of light-activated ChR2 currents, as validated by patch-clamp experiments (see Supplemental Fig. S1). Each pulse was repeated 60 times at a very low rate (0.2 Hz), thus reducing the likelihood of short-and long-term plastic changes in neuronal responses 43 .
Regardless of the pulse duration, the evoked network response outlasted the stimulus and exhibited an early and a late component, as apparent in the shape of the peri-stimulus histogram (PSTH) of the spike times ( Fig. 2A). Control cultures displayed as expected no evoked response (not shown). Similar to electrically evoked responses 37 , the two components were related to the direct neuronal activation and the reverberating interactions mediated by recurrent connectivity. When the intensity and not the duration of the light pulse were changed, considerably less stable and reproducible late responses were observed (not shown), correlating with weaker early responses.
The direct component, occurring in the first few milliseconds after the stimulus onset, originated from the intense sudden depolarization caused by the ChR2 activation. It had low temporal jitter across microelectrodes and repetitions, and was followed by a transient decrease of activity. The reverberating response arose several tens of milliseconds after (i.e. 30-50 ms), as a renewed progressive build up and then exhaustion of neuronal firing. It lasted up to several hundreds of milliseconds and reflected in part the slow deactivation kinetics of ChR2 LC-TC, and in part the effect of recurrent network interactions. As a consequence of the high reliability of the first evoked 1-2 spikes/electrode, observed intracellularly in individual neurons (not shown), the early peak in the PSTH was rather sharp and highly reproducible, with its amplitude and latency not significantly affected by the total pharmacological blockade of synaptic transmission. Amplitude and latency of the early peak in the PSTH did not significantly correlate with the light pulse duration (Fig. 2B). Instead, the peak amplitude and latency of the reverberating response were affected by the total pharmacological blockade of synaptic transmission (not shown) and significantly correlated with the light pulse duration (Fig. 2C,D; r = 0.42 and r = − 0.36, respectively, p < 0.0001): the longer the pulse, the weaker and more delayed the reverberating response. Differences in peak amplitudes and latencies, upon varying pulse duration, were highly significant (p < 0.01) across the entire dataset (n = 39 MEAs), but only for the reverberating responses, not for the early responses.
Examining the evoked responses in more detail, employing higher temporal resolution (i.e. 1 ms bins), we also found signatures of global oscillatory activity in the PSTHs, with striking similarity to physiological gamma-range oscillations. Figure 3A illustrates, for a representative MEA, the spikes evoked by the stimulus as a raster plot of their times of occurrence and its corresponding PSTH, across 59 microelectrodes and 60 repetitions. The estimate of the time-varying power spectrum of the PSTH (Fig. 3B) clearly revealed the appearance of a fading oscillation with frequencies in the gamma range (i.e. ~[40; 100] cycles/sec), as conventionally measured 75 ms after the stimulus onset, and qualified as dominant in terms of signal-to-noise ratio (see the Methods). Across all MEAs tested, oscillations completely faded away ~200 ms after the onset of the stimulus, and generally slowed down at a rate of ~0.3cycles/sec per msec over time. Most importantly, the frequency of the oscillations significantly correlated with the pulse duration (r = 0.55, p < 0.0001): the longer the pulse, the faster the dominant oscillation (Fig. 3C,D).
In addition, an initially low variability of the spike count, estimated by the Fano factor at each microelectrode across repetitions 44 , immediately raised to higher values (i.e., from 0.4 to 1) within 20-50 ms from the stimulus onset (not shown). Moreover, the occurrence time of the spikes recorded at individual microelectrodes showed no entrainment with the oscillation cycles of the PSTH over time, but rather irregular discharges with spike rates generally below 50 spikes/sec (not shown). Overall, these results suggest that excitatory-inhibitory feedback, rather than a beating-wave arising from spike synchronization drift over time, is the prevailing mechanism of the observed oscillations.
Similarity of evoked and spontaneous network events. Inspired by a previous report 45 , we analysed in more detail the time course of the population mean firing rate during spontaneous network-wide events (i.e. network bursts - Fig. 1B) 34,42 , in the absence of any light stimulation. While for evoked responses the spike times were related to the stimulus onset and their histograms averaged accordingly across repetitions ( Fig. 2A), for spontaneous events no temporal reference frame exists by definition. We thus rigidly shifted all spikes, detected simultaneously across the microelectrodes during each spontaneous burst, by a common interval adapted to maximize the firing rate time course similarity across bursts 45 (see the Methods). Unexpected temporal features then emerged in the average firing rate, which could not be revealed by routine analysis methods based on firing rate peaks 42 or threshold-crossing alignment 42 of spontaneous bursts. We found that some but not all MEAs showed prominent signatures of spontaneous transient oscillations in the time course of the firing rate within a network burst, occurring in the same physiological gamma frequency of evoked rhythms (Fig. 4), similarly to the light evoked responses. GABA A receptors are necessary for the evoked global gamma oscillations. In order to investigate the synaptic mechanisms of the evoked oscillations, we pharmacologically abolished fast inhibitory synaptic transmission. As inhibition is necessary for gamma-range rhythmogenesis in vivo 39 , we tested the hypothesis that excitation-inhibition interplay sustains light-evoked oscillations in vitro. When a competitive selective blocker of GABA A postsynaptic receptors was bath-applied, not only the spontaneous network activity significantly increased in terms of firing and burst rate (Fig. 5A) as described 34,46,47 . In addition, while light-evoked responses still consisted of an early and a late component, no oscillatory activity occurred as no single peak in the power spectrum of the PSTHs qualified as a dominant frequency (n = 17 MEAs) (Fig. 5B,C). Hence, the pharmacological manipulation of network activity corroborates the interplay between excitatory and inhibitory neurons, necessary for the light-evoked network oscillations.
Oscillations in silico by the excitatory-inhibitory feedback loop. In order to further investigate the consequences of the excitatory-inhibitory feedback loop, we considered a classic firing-rate mathematical  (Fig. 1A). In all experiments, the time course of the PSTH was stereotyped and consisted of an early and late phase, similarly to what was described for extracellular electrical stimulation. While intensity was fixed, the duration of light pulses was changed systematically and it was unexpectedly found to significantly modulate the late (C,D) but not the early phase of the PSTH (B). Longer pulses significantly (p < 0.001) delayed (C) and weakened (D) the late response peak, without affecting the early phase (circles in A). Evaluating the non-parametric Kendall's rank coefficient, a significant correlation (p < 0.0001) was only found between pulse duration and late-peak latency (0.47) or between pulse duration and late-peak amplitude (− 0.36). Panels B-D summarize the results over 39 MEAs, while the insets show averages across several MEAs, over five distinct sets of sister cultures, remarking the consistency of the findings across different biological preparations.  model 48 . The in vitro cortical networks considered here are composed of both excitatory and inhibitory neurons, showing a prominent pattern of recurrent synaptic connectivity 34 . We therefore defined a model composed of an excitatory population and an inhibitory population with recurrent excitation and feedback inhibition (Fig. 6A), and we analysed and numerically simulated its response to an external constant input. As excitatory neurons express ChR2 in the experiments, only the excitatory population was activated by the external input. After . This range of values for J 0 exists if the feedback inhibition is sufficiently strong, i.e. α > − J ei i 1 (Fig. 6C). Under these hypotheses, oscillations occur at a frequency f, given as the imaginary part of the conjugate complex eigenvalues of the connectivity matrix A (eq. 8), irrespective of the amplitude and duration of the external constant input 49 : implies that feedback inhibition is necessary in our simplified in silico scenario (Fig. 6C), as for J ei = 0 the term below the square root would be negative, and supports the interpretation of the experimental results obtained under blockade of GABA A . Importantly, any (unmodelled) mechanism that should increase the efficacy of excitatory synapses J 0 , e.g. as in short-term synaptic facilitation, would be sufficient to increase f and thus speed up the global oscillations (Fig. 6B). Similarly, any (unmodelled) mechanism that should decrease J 0 , e.g. as a recovery to its resting value following an intense facilitation, would account for the progressive slowdown of the oscillations through time.
Short-term facilitation of synaptic release probability at excitatory synapses. The insights obtained with the model suggest that light, and its duration, might proportionally affect the excitatory synaptic efficacy. We therefore hypothesized that the link between light duration and oscillation frequency might involve the influx and transient accumulation of Ca 2+ in the excitatory synaptic boutons, as in the residual Ca 2+ hypothesis for short-term synaptic facilitation 50,51 . This hypothetical link, occurring by the activation of ChR2 LC-TC and voltage-gated calcium-channels, would lead to a (transient) increase of the action potential-dependent excitatory synaptic efficacy (i.e. analogous to J 0 in eq. 1) and it would compete with the endogenous Ca 2+ recovery processes.
That the synaptic release probability should have been reversibly affected by the light stimulus, and depend linearly on its duration was a prediction of our model.
Intending to test our prediction, we performed whole-cell intracellular recordings from the soma of neurons not expressing ChR2 LC-TC (n = 5 cells), by selecting mCherry-negative cells through videomicroscopy. Then, Figure 5. Inhibition is necessary for light-evoked oscillations. The pharmacological blockade of GABA A receptors disinhibited the network and altered quantitatively but not qualitatively the spontaneous activity, significantly (p < 0.05) increasing the burst rate and decreasing the burst duration (A). However, light-evoked responses were qualitatively affected by network disinhibition as they did not contain oscillatory components (n = 17 MEAs), as revealed by the power spectrum analysis (see Fig. 3) and the lack of a dominant oscillation frequency (B,C).
Scientific RepoRts | 6:24701 | DOI: 10.1038/srep24701 Figure 6. A firing rate model links evoked oscillations to excitatory-inhibitory interactions. By a firing-rate mathematical model (A) we qualitatively reproduced oscillatory transient activity in silico (B). Recruited by the external stimulation of excitatory neurons, the negative feedback of the inhibitory loop generated a well-known alternating pull able to decrease neuronal firing. The red and green traces represent qualitatively the mean firing rate of inhibitory and excitatory neurons, respectively. Altering the efficacy of synaptic connections in a series of numerical simulations (B), resulted in different frequencies of the emerging network rhythm: increasing the strengths of excitatory connections (i.e. J ie and J ee ) increased the frequency of the oscillations. As in the experiments (Fig. 5), inhibition was necessary for the oscillations, as they completely vanished upon removal of the inhibitory feedback (C, J ei = 0). However, the model could account per se neither for the dependency of the oscillation frequency on the duration of the light pulse, nor for its drift over time. As a small increase in the average excitatory synaptic efficacy is sufficient in silico to speed up oscillations, the duration of the light pulse might indirectly affect synaptic efficacy in the experiments, as in short-term facilitation.
Scientific RepoRts | 6:24701 | DOI: 10.1038/srep24701 we blocked action potential (AP) initiation by bath application of a selective blocker of voltage-gated sodium channels (i.e. tetrodotoxin, TTX). Recording in voltage-clamp, we quantified number and instantaneous rate of AP-independent excitatory postsynaptic events 52,53 . Such events are the result of neurotransmitter release, occurring in connected presynaptic boutons spontaneously or reflecting a transient increase in the free Ca 2+ concentration in the boutons 54,55 . We found that immediately after the light stimulus, the presynaptic release probability abruptly and reversibly increased (Fig. 7C), and it was linearly correlated at its peak with the duration of the light pulse (Pearson's r = 0.7, p < 0.0001). This observation is compatible with the residual Ca 2+ accumulation in the synaptic terminal. In addition, the presynaptic release probability reversed significantly over time, exponentially decaying with a time constant of ~200 ms (Fig. 7D), similarly to the recovery from short-term facilitation.

Discussion
Our results show that dissociated neuronal cultures are an interesting model for the study of gamma-range oscillations and are useful to investigate some non-trivial consequences of wide-field stimulation in recurrently connected neurons. Through this experimental model, we report for the first time how to robustly evoke in vitro an oscillatory reverberating network response, and reproducibly modulate its frequency, within a range of physiological importance. Many earlier studies reported on spontaneous and evoked oscillations in vitro 45,56-58 , but none described a fine experimental control of its features. This might prove very relevant for in vitro models of (dys)functional cortical microcircuitry 59,60 and for in vivo optogenetic manipulations of cortical cell assemblies.
We interpreted the emergence of gamma-range damped oscillations from the known interplay between excitation and inhibition, supported by pharmacological experiments and mathematical modeling. An alternative interpretation is offered from computational studies: the firing rate of a population of weakly-coupled neurons shows damped oscillations when all cells respond simultaneously and regularly to the same strong stimulus onset [61][62][63] . Roughly the opposite of a beating-wave, the damped oscillations observed in our experiments might be the consequence of jitter and progressive random shifts of otherwise initially synchronous and regular spike trains. The onset of the light pulse might transiently reset and de facto synchronize the firing of the majority of the cells. Neurons would then start to fire regularly and roughly synchronously, while progressively going out of phase due to distinct membrane properties, heterogeneous ChR2 expression, synaptic interactions, etc.
In such a scenario however, single-neuron firing rates would roughly match the global rhythm of PSTH, with an action potential per cycle. Then, the blockade of feedback synaptic inhibition would not necessarily disrupt the oscillations. For these reasons, and for the physiological range of rhythms observed, the theory of sparsely synchronized oscillations proposed in excitatory-inhibitory network architectures to explain the emergence of gamma oscillations in vivo 39,64,65 , might better describe our in vitro data. In fact, the pharmacological blockade of GABA A receptors completely abolished the oscillations in the PSTH, the spike count across trials showed high variability, and the temporal pattern of firing detected at single microelectrodes appeared to be not completely entrained with the oscillation cycles.
We note however that the blockade of GABA A receptors also removes an inhibitory tone in the excitatory neurons, steepening their frequency-current curves and depolarising their resting potential. It is then possible that these two effects might alter the conditions for oscillations emergence, in a scenario where excitation alone generates the rhythm. Future experiments, mimicking the inhibitory tone in disinhibited networks by optogenetics (e.g. expressing ChR2 and Archaerhodopsin in excitatory neurons), complemented by an accurate biophysical spiking neural model, will be needed to explore the validity of this alternative hypothesis.
In summary, the activity we observed in vitro captures many of the features of in vivo gamma-range rhythms, such as its short-life, the requirement for an intact excitatory-inhibitory loop, and concurrence with irregular single-neuron firing 39 .
To further strengthen this possible conclusion, it would be intriguing in the future to measure the dynamical response 66 of both mCherry-positive/negative cells in vitro, as well as of short-term synaptic plasticity 39,64,65 . Combining these quantitative descriptions together, as described by Wang 65 , one could predict whether or not gamma-range global rhythm can be sustained by recurrent connectivity at the observed firing rates. This might serve as an experimental validation of the theory of sparsely synchronized oscillations 39,64,65 , ultimately linking intrinsic and synaptic properties to emerging oscillations in the same in vitro biological preparation.
An important finding of our experiments is the proposed causal relationship between duration of light pulses and excitatory synaptic efficacy. While the probability of synaptic release was already shown to be altered upon optogenetic stimulation 67 , here we describe a network-level correlate of this effect for the first time. All our experimental results are consistent with our hypothesis of light-induced short-term facilitation of excitatory synaptic efficacy. This requires the reasonable assumption that the transient and reversible increase in Ca 2+ , observed in terms of an enhanced quantal synaptic release, affects the action potential-dependent synaptic release probability, as in the residual Ca 2+ hypothesis for explaining activity-dependent short-term facilitation 50,51 . The impact of ChR2 on Ca 2+ transients has already been documented postsynaptically in dendrites and spines 67,68 , although linked to voltage-gated calcium-channels activation. Here we advanced the hypothesis that direct calcium influx through ChR2 at presynaptic boutons is responsible for the observed phenomenon.
Ultimately, a full spiking neuron model, including a modified version of the Tsodyks-Markram model of short-term synaptic plasticity 69 would be relevant to test the ChR2 LC-TC contribution on synaptic variables. This model would extend the very basic predictions offered by our analysis of the mean-field model ( Fig. 6; eq. 1). Nonetheless, the rate model we examined was simple enough to test the excitation-inhibition feedback loop hypothesis for the oscillations, as well as to explicitly suggest a candidate observable for the additional intracellular experiments of Fig. 7A,B.
An alternative explanation of the relation between the light pulse duration and oscillation frequency might involve the transient alteration of the excitability of (excitatory) neurons expressing ChR2. In fact, as ChR2 LC-TC takes time to deactivate, the total effective membrane ionic conductance in those neurons is larger immediately after the stimulation than before it. This effect would however not account for the increase in the frequency with an increase in the duration of the light pulses. An increase in the effective membrane conductance is known to decrease the slope of single-cell input-output transfer function 70 , thus leading to opposite consequences than those observed, for the frequency f of oscillation (see in eq. 1). In addition, while our hypothesis of presynaptic Ca 2+ accumulation explains why progressively longer pulses result in faster oscillations, a stimulus-dependent change in the total membrane conductance per se would not explain any temporal integration of the stimulus Scientific RepoRts | 6:24701 | DOI: 10.1038/srep24701 feature. Its value would be dependent only on the time since the onset of ChR2 activation, irrespective of the pulse duration.
We also remark that the observed slowdown of the oscillation instantaneous frequency (Fig. 3B) is consistent with our hypothesis of a transient accumulation of Ca 2+ , and its expected depletion during a subsequent recovery phase, as well as with the short-term synaptic depression. As mCherry was expressed over the entire neuronal morphology (Fig. 1D), the spatially unrestricted expression of ChR2, presumably even at excitatory synaptic boutons, supports the Ca 2+ -accumulation hypothesis. Thus, as the brief light pulse instantaneously primed a synapse, subsequent (endogenous) presynaptic action potentials might find easier to release from the same synapse, as in short-term (spiking) activity-dependent facilitation 50,51 .
Our intracellular recordings performed under TTX (Fig. 7B,C) however represent only indirect experimental evidence of the presynaptic facilitation of synaptic efficacy. Intracellular recordings, performed simultaneously from a pair of synaptically connected neurons, would be necessary to quantify whether action potential-dependent synaptic release appears potentiated immediately after delivering light stimulation. While wide-field photostimulation might not be entirely appropriate in that context, a similar experiment might conclusively prove that ChR2 affects recurrent synaptic transmission by altering synaptic efficacy. In addition, with recently developed opsins (e.g. Chronos 71 ), with altered Ca 2+ permeability or faster closing kinetics, we predict that the dependency of the oscillation frequency on the light stimulation duration should change dramatically.
We finally note that the validity of our linear analysis of the mathematical model does not extend to very low firing rate regimes, due to the non-linear response properties of neurons near their firing threshold. In such a regime, oscillations take a more complicated form and their frequency is also determined by the balance between excitation and inhibition.
The finding of a signature of gamma rhythms also during spontaneous bursting, in some but not all MEAs (Fig. 4), is also of great significance to interpret the light activated responses. Despite based on very short (i.e. 5 mins) recordings for each MEA, originally meant as a viability control, this finding generalizes to mature large-scale networks what was first reported in developing small neuronal circuits 45 : an "innate" mode of operation of recurrent neuronal circuits. In addition, our observation indirectly confirms a broad heterogeneity in the excitation-inhibition ratio in cortical networks, reflected in the variability of the spontaneous rhythms' frequency from ~50 to ~100 cycles/sec (Fig. 4) (or the lack of it) 72 . In addition, transitions were sometimes observed (e.g. Fig. 4A, lower panel) between non-oscillating and oscillating modes of intra-burst activity, consistent with the broad range of time-scales for single-neuron and network excitability 72 . Finally, our results further support the dynamical character of such an excitation-inhibition ratio 72 , as it can be readily manipulated by internal or external inputs (i.e. light pulses).
Overall, in line with many previous studies 4,8,14,73 , our results stress the necessity to take into account the possible impact on recurrent connectivity when designing optogenetic experimental protocols. Moreover, we emphasize the importance of further developments in the field, particularly confining the opsin expression in specific regions of a neuron 74 , based on exact experimental needs. In fact, the light-induced short-term facilitation of excitatory synaptic efficacy might prove to be a bug for some applications and not a feature. Ultimately, besides further investigating and clarifying the impact of optogenetic manipulations in large-scale recurrent cortical networks, particularly those affecting synaptic release probability, our work provides a means to manipulate the network responses, for the benefit of future principled design of closed-loop control paradigms.
Large-scale cultured networks are acknowledged as an experimental preparation, where structural, dynamical, and plastic properties of in vivo cortical networks may be effectively investigated 22,32,34,38,42,[75][76][77][78][79][80][81] . Due to their inherent 2-dimensional geometry, lower effective cell packing density, and random arrangement of cells' morphology in the plane, they cannot capture faithfully the features of electrical local field potentials (LFPs) as recorded in 3-dimensional layered structures as brain slices and in vivo 82 . While our (evoked or spontaneous) in vitro gamma oscillations occur prominently in the network-wide probability of firing (Fig. 3A), intact brain rhythms are often reflected also in the spatially-integrated activity of LFPs. As we examined the low frequency components of single-trial MEA raw recordings, a weak signature of the global firing rhythm could be indeed observed (Supplemental Fig. S2) although with orders of magnitude worse signal-to-noise ratio and more difficult interpretation 83 than the spike PSTHs.
In conclusion, the emergence of oscillations at physiological gamma frequency and their modulation by the selective activation of excitatory neurons makes large-scale cultured networks, MEAs, and optogenetics relevant for studying and manipulating network phenomena under highly-controlled conditions (see also [84][85][86], as well as for validating theoretical models of population dynamics with simplified hypotheses on unstructured connection topology.

Materials and Methods
Viral vectors production. To replace the cytomegalovirus (CMV) promoter in the Adeno-Associated Viral vector (AAV) transfer plasmid 66 , a 364 bp fragment of the mouse α -subunit of Ca 2+ / Calmodulin-dependent protein kinase II promoter (CaMKIIα ) was PCR-amplified from plasmid 20944 (Addgene, UK), using Fw-AAAGCTAGCACTTGTGGACTAAGTTTGTTCACATCCC-and Rev-AAAAGCGCTGATATCGCTGCCCCCAGAACTAGGGGCCACTCG-as primers. This PCR fragment was ligated in the AAV transfer plasmid, using NheI and Eco47 III, after the CMV promoter was cut out. eGFP was then replaced by ChR2 L132C/T159C-mCherry. The opsin coding sequence was cloned by PCR amplification and the final transfer plasmid was sequence-verified prior to viral vectors production. AAV2/7 viral vectors were produced at the Leuven Viral Vector Core, as described earlier 66 . Briefly, HEK 293T cells (ATCC, Manassas, VA, USA) were seeded in the Hyperflask at 1.2 × 10 8 cells per viral vector production, in Dulbecco's Modified Eagle's Medium (DMEM) with 2% Fetal Calf Serum. The following day the medium was replaced by serum-free Optimem and cells were transfected with the pAAV transfer plasmid, the rep/cap plasmid, and the pAd.DELTA.
Scientific RepoRts | 6:24701 | DOI: 10.1038/srep24701 F6 plasmid in a 1:1:1 ratio. The supernatant was harvested five days after transient transfection, concentrated using tangential flow filtration, and purified using an iodixanol step gradient. Aliquots were stored at − 80 °C and quantified, using real-time PCR, as Genome Copies per ml (GC/ml): titers varied between 9.8 × 10 11 GC/ml and 1.6 × 10 12 GC/ml. Neuronal culturing and transduction. Primary cultures of mammalian neurons, dissociated from the postnatal rat neocortex, were prepared as in 30 , in accordance with international and institutional guidelines on animal welfare. All experimental protocols were specifically approved by the Ethical Committee of the Department of Biomedical Sciences of the University of Antwerp (permission n. 2011_87), and licensed by the Belgian Animal, Plant and Food Directorate-General of the Federal Department of Public Health, Safety of the Food Chain and the Environment (license n. LA1100469). Briefly, cerebral cortices (excluding the hippocampus) were removed from the brains of Wistar pups (P0), flushed with ice-cold phosphate buffered saline (PBS) containing 20 mM glucose, and roughly minced with a blade. Hemi-cortices were enzymatically digested by adding 0.05 mg/ml trypsin to the solution and by gentle agitation, for 15 min at 37 °C (GFL 1083, Gesellschaft für Labortechnik mbH, Burgwedel, Germany). After stopping the digestion by adding 2 ml of heat-inactivated Normal Horse Serum (NHS), the tissue fragments were left to sediment at room temperature and then mechanically triturated by a 10 ml Falcon pipette, filtered through a nylon monofilament cell strainer (40 μm, #352340, BD Falcon, Franklin Lakes, NJ, USA), and centrifuged at room-temperature (220 g for 5 min; GS-6, Beckmann-Coulter, Brea, CA, USA). The resulting cell pellet was resuspended in culture medium, composed of 4 ml Minimum Essential Medium (MEM), 5% NHS, 50 μg/ml gentamycin, 0.1 mM L-glutamine, and 20 mM sucrose, and diluted to reach a final surface plating density of 6,500 cell/mm 2 on glass coverslips (15 mm, VD1-0015-Y2MA, Laborimpex, Forest/Vorst, Belgium) and commercial substrate-integrated arrays of microelectrodes (MEAs; Multi Channel Systems GmbH, Reutlingen, Germany).
MEAs with a regular 8 × 8 arrangement of 60 Titanium Nitrate (TiN) microelectrodes, each with 30 μm diameter and 200 μm spacing (60MEA200/30iR-ITO-gr, Multi Channel Systems), were employed. Prior to cell seeding, MEAs and coverslips were coated overnight with polyethyleneimine (0.1% wt/vol in milliQ water at room temperature), and then extensively washed with milliQ water and let air-drying. Seeded MEAs and coverslips were maintained in an incubator for up to 6 weeks at 37 °C in 5% CO 2 and 100% R.H. (5215, Shellab, Cornelius, OR, USA), with MEAs sealed by fluorinated Teflon membranes (Ala-MEA-Mem, Ala Science, Farmingdale, NY, USA) and coverslips stored in Petri dishes.
The culture transduction was performed in vitro five days (DIV5) after seeding, by partly replacing the medium of each MEA or coverslip with a 1:50 dilution of viral particles AAV-CaMKIIα -ChR2(LC-TC)-mCherry in fresh medium. At DIV8, fresh and pre-warmed medium was added to reach 1 ml final volume in MEAs and coverslips. From DIV10 on, half of their culture medium volume was replaced every 2 days with fresh, pre-warmed medium. All reagents were obtained from Sigma-Aldrich (St. Louis, MO, USA) or Life Technologies (Gent, Belgium).
MEA extracellular electrophysiological recordings. Before each experiment, MEAs were transferred from the incubator to an electronic amplifier (MEA-1060-Up-BC, Multi Channel Systems) with 1-3000 Hz bandwidth and amplification factor of 1200, placed inside an electronic-friendly incubator (37 °C, 5% CO 2 , ~20 R.H.). MEAs were left in the incubator to accommodate for 5 min before starting the recordings, and were used to detect non-invasively the spontaneous and (light) evoked neuronal electrical activity. An MCCard A/D board and the MCRack software (Multi Channel Systems) were used to sample (25 kHz/channel), digitize (16 bits), and store on disk the raw voltage traces individually detected by each of the 60 microelectrodes for subsequent analysis. Recordings took place after 28 DIVs, as continuous as well as triggered acquisitions of raw traces. Assessing the health of each MEA, five minutes of spontaneous activity was routinely recorded, before delivering the photo-stimulations (Fig. 1B,C). Light-evoked activity was then recorded across the entire MEA, by monitoring 200 ms proceeding and 1000 ms following the photo-stimulus onset. Intracellular electrophysiological recordings. Glass coverslips were employed for patch-clamp experiments in sister cultures of at least 20 DIV, using an Axon Multiclamp 700B amplifier (Molecular Devices, USA) controlled by the LCG software 67 . Patch pipettes were pulled on a horizontal puller (P97, Sutter, Novato, USA) from filamented borosilicate glass capillaries (World Precision Instrument Inc. (WPI), USA), and had a resistance of 4-5 MΩ, when filled with an intracellular solution containing (in mM): 115 K-gluconate, 20 KCl, 10 HEPES, 4 Mg-ATP, 0.3 Na 2 -GTP, 10 Na 2 -phosphocreatine (pH adjusted to 7.4 with NaOH). An extracellular solution, containing (in mM): 145 NaCl, 4 KCl, 1 MgCl 2 , 2 CaCl 2 , 5 HEPES, 2 Na-pyruvate, 5 glucose (pH adjusted to 7.4 with NaOH), was constantly perfused at a rate of 1 ml/min and at a temperature of 37± 1 °C for each experiment.
Miniature post-synaptic currents (mPSCs) were recorded under voltage-clamp at a hyperpolarized holding potential of − 70 mV, under the whole-cell configuration from the soma of mCherry-negative neurons, after low-pass filtering at 6 kHz and sampling at 30 kHz the acquired traces. Photo-activated responses were also recorded monitoring the membrane currents, while patching mCherry-positive neurons. Repeated stimulation trials, acquiring 2 s preceding and 2 s following each photo-stimulation were performed, as in MEA experiments. Cells with resting potential above − 55 mV, or with a series resistance larger than 25 MΩ, were not included for analysis and discarded. In all the experiments, the Active Electrode Compensation 87 was employed by default through LCG.
Pharmacology. In some MEA experiments, Gabazine (GBZ, SR95531) was bath-applied prior to the electrophysiological recordings at the final concentrations of 30 μM, as a selective competitive blocker of GABA A receptors (Fig. 5). In all intracellular experiments (Fig. 7A,B), 1 μM tetrodotoxin (TTX) was added to the extracellular solution to selectively block sodium channels and thus suppress the generation of action potentials in all neurons. All drugs were obtained from Abcam (Cambridge, England, UK).
Data Analysis and Statistics. For MEA recordings, extracellular voltage waveforms of light-evoked neuronal responses were analyzed on-the-fly by custom scripts written in MATLAB (The MathWorks, Natik, MA, USA), while the package QSpike Tools 69 was employed for analyzing spontaneous activity and extracting standard observables such as spike times, spike rate, burst rate, and burst duration 22,42 . Raw voltage traces recorded at each microelectrode ( Fig. 1A) were zero-phase band-pass filtered (400-3000 Hz) and then processed by an adaptive peak-detection algorithm, based on voltage-threshold crossing and including an elementary sorting of positive and negative amplitude peaks of the corresponding spike waveforms [88][89][90] . The times of occurrence of extracellular spikes, for each recording channel and each repetition, were stored on disk for subsequent analysis together with their corresponding spike waveform. Spike waveforms with negative peak voltage amplitudes were sorted and included in the analysis. Specifically, their time stamps were used to estimate network-wide evoked instantaneous spike response probability, computed as a Peri-Stimulus spike Time Histogram (PSTH) over 5 ms bins ( Fig. 2A) and whose features were subsequently extracted (Fig. 2B,D).
Spontaneous network bursts were also studied, analyzing the corresponding firing probability over time by computing for each burst a spike time histogram (STH) over 5 ms bins. The overall firing rate profile during a burst, averaged over all the bursts recorded from a MEA, was however not obtained arbitrarily aligning the STHs recorded to the time of their peak 42 . Instead, STHs were first arbitrarily shifted in time and then aligned for averaging, by a series of time intervals that maximized their mutual similarity, through maximizing their correlation coefficients 45 .
Frequency-domain analysis was also performed, estimating the spectrogram 49 of each PSTH (estimated over 1 ms bins, Fig. 3A) corresponding to a different photo-stimulus duration (Fig. 3B). In details, the Fast Fourier Transform with a sliding window (50 ms width, 50% overlap) was employed to extract the time-varying spectrum of frequencies contained in the PSTH 91 , averaging over the 60 trials and normalizing across frequencies by the PSTH spectrum of the 200 ms preceding the photo-stimulus onset. In addition, the power spectrum computed 75 ms after the stimulus onset (Fig. 3B,C) was conventionally employed to determine existence and location of the dominant oscillation frequency, defined arbitrarily as the frequency corresponding to the highest peak of the spectrum, which necessarily exceeded the median value of the spectrum over all frequencies by three times (Fig. 3C).
In intracellular recordings, current traces were low-pass filtered at 1.5 kHz and analysed using custom-written MATLAB scripts to count the number or the instantaneous rate of mPSCs.
Data is presented as mean ± standard error of the mean (SEM), except for Figs 1C and 5A, where box plots show the median (horizontal bar) and the mean values (circle), the 25th and 75th percentiles (box's edges), and the most extreme data points (whiskers). Distributions were compared using the Kolmogorov-Smirnov test and only differences with at least p < 0.05 were considered significant. Correlations were assessed by computing the non-parametric Kendall's rank correlation coefficient and its significance level 85 .
Wide-field photo-stimulation. A blue light-emitting diode (LED) (Rebel, Quadica Development, Canada), powered by an externally dimmable DC driver (LuxDrive, Randolph, USA), was employed to deliver wide-field photo-stimulations. The LED was further equipped with a parabolic reflector lens, to collimate the light beam and distribute the emitted power on the bottom of MEAs and coverslips uniformly (Fig. 1A).
In MEA experiments, photo-stimuli were delivered at full power (i.e. 2 mW/mm 2 , measured at the sample by a calibrated photodiode; 818-ST2-UV, Newport Spectra-Physics, Netherlands) by a voltage-controlled stimulus generator (STG1002, Multi Channel Systems) connected to the DC driver. The generator was programmed over USB to define timing and waveform of each stimulus. 60 identical square pulses were consecutively delivered at a 0.2 Hz repetition rate, and the specific effect of the pulse duration [0.1, 0.5, 1, 2, 5, 10, 20] ms was investigated. Thus, stimulation sessions differed only in the pulse durations, whose value was chosen in a shuffled order across successive sessions.
In the intracellular experiments, an analog output of the acquisition D/A board was connected to the LED DC driver and the LCG software employed to deliver the photo-stimuli as in MEA experiments, repeated 20 times per condition at 0.2 Hz.
In our experiments, the pulse duration T was varied while keeping its intensity L max to full power to improve proportionality between stimulus and net inward charge, as tested by voltage-clamp experiments (Supplemental Fig. S1). In fact, the equation describing the fraction x of open ChR2 channels, from a minimal 2-state (non-inactivating) kinetic biophysical model, predicts a net charge proportional to T, only for large intensities L max .
Firing-rate model. A Wilson and Cowan rate-based mathematical model was defined and numerically simulated 48,92 , to review the consequences of recurrent and feed-back interactions between an excitatory (e) and an inhibitory (i) neuronal population. Each population was assumed to be composed of identical and indistinguishable neurons and thus collectively described by a single state variable (i.e. h e (t) and h i (t)). This variable approximated the total average incoming synaptic inputs to a generic unit of each population. The time-varying output mean firing rate of a population (i.e. E and I) was then assumed, for simplicity, to depend only on its total synaptic input, as in a threshold-linear frequency-current curve (eq. 4): e e e ii i where + [ ] indicates the positive part of its argument, θ a minimal activation threshold, and α the slope of the frequency-current response function. In the case of a network topology including recurrent excitation and reciprocal excitation and inhibition (Fig. 6A), the following coupled (non-linear) differential equations fully describe the dynamics of the system: e e e e e e i 0 , , ee ei ie ii indicate the average synaptic efficacies of recurrent excitation, reciprocal excitation and inhibition, and mutual inhibition, respectively. For simplicity of the subsequent analysis, we assumed no reciprocal inhibition (i.e. = J 0 ii , see Fig. 6A). E 0 represents the time-varying external input current, induced by ChR2 activation, and τ e and τ i are the time constants of decay of the synaptic currents in the lack of presynaptic activity, for each population respectively. Taking into account the rather slow kinetics of ChR2 LC-TC, E t ( ) 0 was approximated as a brief square pulse followed by a much longer decaying exponential function. By a change of variables, θ = − x h e e and θ = − y h i i , and assuming that both h e and h i at a given moment are larger than θ e and θ i , respectively, eqs 5 and 6 can be rewritten as a linear system of differential equations, using linear algebra notation: where ii , determine the eigenvalues of the matrix A and thus the dynamics of the output network response (e.g. x(t)) to the external stimulus (i.e. E 0 (t)), being exponential or oscillatory over time 93 . The parameters employed for the simulations of