Theta rhythm-like bidirectional cycling dynamics of living neuronal networks in vitro

The phenomena of synchronization, rhythmogenesis and coherence observed in brain networks are believed to be a dynamic substrate for cognitive functions such as learning and memory. However, researchers are still debating whether the rhythmic activity emerges from the network morphology that developed during neurogenesis or as a result of neuronal dynamics achieved under certain conditions. In the present study, we observed self-organized spiking activity that converged to long, complex and rhythmically repeated superbursts in neural networks formed by mature hippocampal cultures with a high cellular density. The superburst lasted for tens of seconds and consisted of hundreds of short (50–100 ms) small bursts with a high spiking rate of 139.0 ± 78.6 Hz that is associated with high-frequency oscillations in the hippocampus. In turn, the bursting frequency represents a theta rhythm (11.2 ± 1.5 Hz). The distribution of spikes within the bursts was non-random, representing a set of well-defined spatio-temporal base patterns or motifs. The long superburst was classified into two types. Each type was associated with a unique direction of spike propagation and, hence, was encoded by a binary sequence with random switching between the two “functional” states. The precisely structured bidirectional rhythmic activity that developed in self-organizing cultured networks was quite similar to the activity observed in the in vivo experiments.


Introduction
Synchronization and the interplay between excitation and inhibition in neural networks play crucial roles in the organization of rhythmic activity in the brain [1][2][3][4][5]. Rhythmic oscillatory activity with various frequencies represents a multi-clock substrate for cognitive function, memory and sleep [6,7]. However, researchers still question whether the rhythmicity emerges from the specific network morphology that develops during neurogenesis [7][8][9][10][11] or it is generated spontaneously due to nonlinear network dynamics mediated by an interplay between excitation and inhibition that is sustained by a homeostatic balance [12][13][14]. An answer to this a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 Superbursts in neuronal cultures represent highly coherent spatio-temporal spiking activity patterns that spontaneously develop in originally non-structured networks due to self-organization and plasticity [22].
One of the proposed mechanisms for superbursts is the interaction between the activity of excitatory and inhibitory neurons in culture. In particular, the addition of inhibitory cells from the striatum to the hippocampal cultures was used to study their impact on burst dynamics. An increase in the inhibitory cell fraction in neuronal cultures from 20% to 56% significantly increases the number of small bursts in the superburst structure [47]. The important role of GABAergic neurons in the generation of bursting activity has also been confirmed in network mathematical models in which GABAergic neurons were involved in generating the small bursts in subsequent superburst [46]. Superbursting activity has been increased by treating cultures with blockers of inhibitory synapses (bicuculline and picrotoxin) or a low concentration of Mg 2+ ions [18,48] and has been decreased by treating cells with a combination of Na + channel blockers and picrotoxin [18].
The investigation of spiking patterns in the superbursts has revealed remarkably precise repetition of the internal bursting sequence [22]. Superbursts appear with irregular intervals, but their internal structure contains small bursts with highly regular and reproducible activation patterns that persist for hours or days [26]. In another study of cortical cultures during the mature stage (4)(5)(6) week in vitro), definite motifs observed in the burst activation pattern corresponded to a specific oscillation phase during the ultra-slow oscillations (<0.01 Hz) [12]. These ultra-slow oscillations also represent superbursts with regular intervals. Spiking pattern motifs were found to be strongly conserved across multiple oscillation cycles, repeating themselves with high spatio-temporal precision.
According to one recent mathematical model,high neurite and synapse densities may also influence small bursts in the superburst subsequences [46]. Stable rhythmic activity in the form of propagating synchronized bursts over several minutes was induced in cortical cultures treated with inhibitors of GABAergic synaptic transmission in another study [13]. This periodic synchronized activity on the 3-4 second time scale was only observed at the boundaries of the culture. Thus, the excitatory-inhibitory balance should be an important parameter for generating stable and reproducible synchronized activity in networks.
In the present study, we observed long superbursting activity with well-defined and reproducible temporal dynamics in spontaneously developed hippocampal neurons cultured on a microelectrode array (MEA). Regarding the electrophysiological activity, we observed long (up to 30 seconds) superbursts consisting of subsequences of (up to hundreds of) highly reproducible short bursts in the centre of the cultured network. The spiking frequency in the bursts was 139.0 ± 78.6 Hz, and the interburst interval ranged from 100-150 ms (11.2 ± 1.5 Hz), which resembled unique hippocampal activity under in vivo conditions [49]. Spike propagation pathways during short bursts in all long superbursts were aligned along two major spatial directions. The long superburst was encoded into two types, each associated with a definite orientation of spike propagation. The orientation was switched during each single superburst; in the subsequent superburst, the orientation was determined randomly, with a probability of switching to the next orientation of approximately 50%. Therefore, the superburst time sequence was encoded by binary symbols reflecting the spontaneous activation of two dominant spike propagation patterns selected in mature networks of cultured hippocampal cells. Notably, this well-organized rhythmic activity emerged spontaneously in mature culture networks in vitro without any specific stimulation and afferentation. We believe that these selforganizing dynamics observed during culture development led to a certain excitatory-inhibitory balance, where cycling dynamics serve as a homeostatic functional state of the culture. Moreover, these findings also suggest a possible mechanism by which different functional states (i.e., vortices and synchronized epileptic-like discharges) may spontaneously appear in brain networks in the absence of stimulation.

Cell culture
Hippocampal cells were dissociated from embryonic mice (E18) and plated on microelectrode arrays (MEAs) pre-treated with adhesion promoting molecules of polyethyleneimine (Sigma P3143) with a final density of approximately 15,000-20,000 cells/mm 2 (Fig 1D). Our cultures were composed of 4-5 layers of the cells. The mice used in our study were received from Institute of Bioorganic Chemistry Pushchino, Moscow Region, Russia. C57Bl/6 mice were euthanized via cervical dislocation according to protocols approved by the National Ministry of Public Health for the care and use of laboratory animals. The protocol was approved by the Committee on the Ethics of Animal Experiments of the Nizhny Novgorod State Medical Academy (Permit Number: 9-25.09.2014). All efforts were made to minimize suffering. Embryos were removed and decapitated. The entire hippocampus was dissected under sterile conditions. The cortex, whole medulla and the lower part of the pons were excluded during the dissection. Hippocampi were cut in Ca 2+ -and Mg 2+ -free phosphate-buffered saline (PBS-minus). After enzymatic digestion for 25 min using 0.25% trypsin (Invitrogen 25200-056) at 37˚C, cells were separated by trituration (10 passes) using a 1 ml pipette tip. Next, the solution was https://doi.org/10.1371/journal.pone.0192468.g001 centrifuged at 1500 g for 1.5 min, and the cell pellet was immediately re-suspended in Neurobasal culture medium (Invitrogen 21103-049) with B27 (Invitrogen 17504-044), glutamine (Invitrogen 25030-024) and 10% fetal calf serum (PanEco 055). The dissociated cells were seeded in a 30 μl droplet covering the centre of the culture dish with 1 mm 2 electrode region of the MEA. It resulted in a culture 6-7 mm in diameter. After the cells had adhered (usually in 2 hrs), the dishes were filled with 1 ml Neurobasal medium (NBM) supplemented with B-27 and 0.5 mM glutamine with 10% fetal calf serum. After 24 hrs, the plating medium was replaced with a medium containing NBM 2% B-27 and 1% glutamine and 0.5% fetal calf serum but with no antibiotics or antimycotics. Glial growth was not suppressed, given that glial cells are essential for long-term culture health. Half of the medium was replaced every 2 days. The cells were cultured under constant conditions of 35.5˚C, 5% CO 2 and 95% air at saturating humidity in a cell culture incubator (MCO-18AIC, SANYO).
Phase-contrast images of cultures were taken weekly to record the status of the culture using a Leica DMIL HC (Germany) inverted microscope with 10/0.2 Ph1 objectives. Experiments were conducted when the cultures had been grown for 3-5 weeks in vitro.

Electrophysiological methods
Extracellular potentials were collected using 59 planar TiN electrodes integrated into the USB-MEA-120 system (Multichannel system, Germany). The microelectrode arrays (MEA) had 59 electrodes (8x8 grid) with diameter of 30 μm and spaced 200 μm apart ( Fig 1A). Data were recorded simultaneously from 59 channels at a sampling rate of 20 kHz/channel. All signal analysis and statistics were performed using the custom-made software Meaman in Matlab (Mathworks, USA).

Spike detection
The detection of recorded spikes (Fig 1B) was implemented using threshold calculation: where σ = median (|x| / 0.6745), which was the estimate of the median normalized to standard deviation of a signal with no spikes (see [50] for more details), x is the band-pass-filtered (0.3-8 KHz) signal, and NS is the spike detection coefficient, which was set to 8. The amplitudes of detected spikes were in the range of 20-200 μV. The minimal interspike interval was set to be 1 ms to avoid the overlapping of neighbouring spikes.

Burst detection
The burst detection method was described in detail in our previous paper [31]. Briefly, we estimated the total spiking rate characteristic, TSR(t), as the number of spikes from all electrodes within each 5 ms time bin. The rapid appearance of a large number of spikes over all electrodes in a small (2 ms) time bin was used as the criterion for burst appearance. Threshold detection was applied to estimate the initiation and termination of the bursts. The burst threshold was set to T Burst = 0.2 × σ TSR , where σ TSR is the standard deviation of TSR(t).
The initiation time of the burst was defined as the burst start time, where TSR(t) was greater than the threshold. Next, the initiation time was adjusted to the first spike from all electrodes after a supra-threshold time. Finally, the time point at which TSR crossed the threshold after the burst start time was defined as the burst end time.
Interburst peak intervals (IBPIs) were calculated as the time interval between adjacent peaks of the spiking rate of detected small bursts. The peaks were separately estimated from the total spike rate (TSR) diagram of each burst. Examples of the peaks and the interval between two bursts are presented in Fig 2B. Such measure of interval between "up" states of the bursts is similar to a measure of oscillation period and frequency in rhythmic activity analysis of in vivo recordings. Each IBPI corresponded to the instantaneous frequency (IF) of each pair of the small bursts. All burst pairs were analysed, and IF datasets were used for further statistical analysis.

Superburst detection
Superbursts and long superbursts in the electrical activity were detected using a previously described method [51]. First, we defined a Gaussian function with an effective width equal to 50 s. Next, that function was iteratively moved from the beginning of the recording to the end using a 10 ms time step, while the cross-correlation of the function with the TSR was calculated for each step. The resulting cross-correlation indicated the amount of synchronized activity (bursts) that was recorded in each 10 s window. We applied a threshold detection algorithm in which the threshold was estimated as the superburst detection accuracy coefficient multiplied by the standard deviation of the calculated cross-correlation to detect superbursts in the spiking activity. The superburst detection accuracy coefficient was estimated empirically and was equal to 0.4. All time points that crossed the threshold were defined as the beginnings and the endings of the superbursts [51].

Burst classification
Superbursts consisted of initiation bursts lasting for 50-100 ms and short small bursts lasting for 30-50 ms. The total number of spikes within initiation bursts ranged from 1000-3000 spikes, whereas each small burst contained 10-500 spikes ( Fig 1F). These two types of bursts were identified using a K-means clustering algorithm.
We analysed activation patterns consisting of first spike timings of the bursts to represent the spatio-temporal properties of all patterns within small bursts. The first spike timing was averaged for each electrode and each small burst. Then, the values from all 60 electrodes were colour coded and plotted on an image using cubic convolution interpolation ( Fig 3F, 3G and 3H). This image represented a gradient map of the burst activation profile.
We introduced a vector field map of activation timings in the culture to represent the spatial properties of spike propagation during burst activation. For each MEA electrode, we calculated a vector whose direction represented the activation-timing gradient around an area of 3 electrodes ( Fig 3F). Notably, the resulting vector field resembled a colour-coded activation pattern. We defined this spatial representation of the activation pattern as a dynamic pattern.
We applied the K-means clustering method to identify motifs of activation patterns in small bursts of all long superbursts. The activation patterns for each small burst consisted of the first spike timings for each of the 59 electrodes of the MEA. This method required a number of clusters to be estimated. First, we estimated two clusters and evaluated cluster separation by calculating the Davies-Bouldin (DB) index [32]. The DB index estimates the ratio between the internal cluster distance and the distance between clusters. Then, the clustering procedure was repeated for various numbers of clusters (2, 3. . .30), and the DB index was estimated ( Fig  4C). The minimum value for the DB index among all tested cluster numbers was 2 or 3, indicating that the activation patterns were optimally clustered into 2 or 3 motifs of small bursts. Small values for the DB index corresponded to compact clusters whose centres were located far from each other. DB values ranging from 0 to 1 indicated "robust" clustering. We also tested the motif separation using the expectation-maximization algorithm (EM clustering). We plotted two principal component coefficients for each activation pattern and highlighted the estimated clusters in different colours to represent the evaluated clusters ( Fig 4D). As in the previous case, we used the DB index to evaluate the optimal number of clusters with EM clustering (Fig 4F). This method was applied to 3 principal component coefficients and divided data into the clusters more accurately (Fig 4H, 4I, 4J and 4K), which were used in subsequent analyses.

Results
First, we analysed the spontaneous activity of the hippocampal cultures. We obtained complex bursting patterns similar to those reported previously in cortical cultures [26]. An example of (D) Distribution of spiking rate frequency from the electrodes with small bursts in one culture. The mean spiking rate was 178.5 Hz. (E) Distribution of IFs (n = 6 cultures). The median burst frequency was 11.2 ± 1.5 Hz (mean ± s.d., n = 6 cultures). (F) Distribution of the spiking rate frequency from the electrodes with small bursts (n = 6 cultures). The spiking rate was 139.0 ± 78.6 Hz (mean±s.d.).
https://doi.org/10.1371/journal.pone.0192468.g002 Examples of the bursts from 3 motifs (C, D, and E) and spatial representation of dynamic patterns (F, G, and H) are shown. Motif #1 was observed in 10.6% of small bursts, motif #2 was observed in 56.4%, and motif #3 was observed in 32.9% of all small bursts in the superburst. Examples of two types of long superbursts (I and J) that were composed of the 3 motifs. The left long superburst consisted of motif #1 bursts (blue markers) and motif #2 bursts (green markers) (F). The second type of long superburst consisted of bursts of motif #3 (red markers, J). the spikes recorded from a single electrode within a small burst is shown in Fig 1C. After 3-4 weeks of culture in vitro, we obtained the activity described as a superburst (Fig 1E). A typical superburst consisted of a sequence of 3-20 small bursts of 50-100 ms in duration and a 50-150 ms interburst interval. During the period of 30-40 DIV, the cultures generated long superbursts that were similar to regular superbursts, but that lasted for 10-30 seconds and consisted of hundreds of regular small bursts. In summary, we analysed 11 cultures from 3 plating experiments and observed long superbursts in 8 cultures. Six cultures generated more than 6 long superbursts for at least 20 minutes, which were included in the statistical analysis. In other cultures, we observed no more than 2 long superbursts. The signals from a single electrode during the long superburst on timescales of 15 and 2 seconds are illustrated in Fig 1A  and 1B. Raster plots of the spiking activity recorded from all 59 electrodes during the superburst and long superburst are shown in Fig 1G. Each point on the raster plot represents the time at which a spike occurred at a particular electrode. The long superbursts were composed of relatively long initiation bursts (50-100 ms) followed by shorter bursts, i.e., the small bursts. The initiation bursts and the small bursts were easily identified by K-means clustering (Fig 1F) using burst firing rate features (see the Methods).
Next, we estimated the frequencies of the small bursts; the initiating bursts and intervals between the long superbursts were excluded from analysis. We estimated instantaneous frequencies of the detected bursts using the interburst peak interval measure (IBPI) (see the Methods) (Fig 2A). Each IBPI corresponded to the instantaneous frequency (IF) of each pair of the small bursts. Notably, the IFs were not normally distributed (Kolmogorov-Smirnov test, p<0.01), and the median value of the IF was estimated for each culture. A typical example of the IF distribution for a small burst sequence in one culture is shown in Fig 2C. Moreover, a median value of the IFs represented a bursting frequency, which was stable and equal to 9.8 Hz. Less than 5% of the IFs ranged from15-30 Hz in presented example. Then, the bursting frequencies were averaged for all cultures, and a mean value was equal to 11.2 ± 1.5 Hz (mean ± standard deviation, n = 6 cultures). An average histogram of the IFs for all 6 cultures is illustrated in Fig 2E demonstrating high reproducibility in various preparations. Most of the IFs were concentrated in the range from 8 to 15 Hz.
We also estimated the mean spiking frequency for each electrode during small bursts exclusively during the intra-burst periods ( Fig 2B, the red rectangle marks the period in a sample burst). For the raster plot presented in Fig 1A, the spiking frequency at most electrodes ranged from 100-300 Hz, and the mean frequency was 178.5 Hz (Fig 2D). On average, the spiking frequency per electrode was 139.0 ± 78.6 Hz (mean ± standard deviation, n = 6 cultures, 264 active electrodes) ( Fig 2F).
Next, we analysed spiking patterns in sequences of small bursts within a long superburst recorded for 30 minutes. The set of the first spikes in the burst recorded from all electrodes was considered the activation pattern [31]. We applied the EM clustering algorithm for 3 principal component features to investigate different motifs (clusters) of activation patterns. This method required the estimation of a number of clusters. First, using two clusters, we evaluated the cluster separation by calculating the DB index (see the Methods). Then, the clustering procedure was repeated for various numbers of clusters (2, 3. . .30), and the DB index was estimated ( Fig 3A). The minimum value for the DB index among all tested numbers of clusters was 3, indicating that the activation patterns were optimally clustered into 3 motifs in the presented raster plot. The activation patterns and profiles of the spiking patterns in the bursts from separate clusters (motifs) represented different sequences of spike occurrence, i.e., different dynamics of the spike propagation (Fig 3B). Note that the difference in the first spike time sequence in the pattern can be visually observed in Fig 3. Representative raster plots for all small bursts from all 3 motifs are shown in Fig 3C, 3D and 3E. We averaged the activation patterns and calculated the dynamic patterns to investigate the spatio-temporal properties of all patterns within each motif (see the Methods ( Fig 3F, 3G and 3H)). The first spike timings from all 59 electrodes were colour coded and plotted on the image created using cubic convolution interpolation. Surprisingly, burst activation was implemented in the form of wave-like spike propagation dynamics with a wide wave front. Arrows represent the gradient of activation times, i.e., the mean direction of spike propagation during burst initiation across each electrode. Eventually, the patterns were organized into a uniform direction in space. Remarkably, motifs #1 and #2 presented similar directions of the activation pattern, from the upper electrodes to the bottom of the MEA, whereas motif #3 presented the opposite direction.
Next, we analysed the sequence in which the motifs appeared in the structure of the long superburst. Some of the long superbursts (Fig 3I) comprised bursts of motifs #1 and #2, whereas other long superbursts in the same recording (Fig 3J) mainly comprised bursts of motif #3.
Notably, the firing rate in the burst sequence (Fig 3J, TSR-total spiking rate) was quite variable, but the sequence of the first spike timings, i.e., the activation patterns, of the bursts remained largely unchanged.
Furthermore, we verified the clustering results for activation patterns using EM clustering for two principal components (PCs) and three PCs and K-means clustering (see the Methods).
We applied all three methods to one data set; we varied the number of clusters and estimated clustering using the DB index as shown in Fig 3A. Using K-means and EM clustering with 2 PCs, the DB index had a minimum at 2 clusters (Fig 4C and 4F), whereas the DB index calculated using EM with 3 PCs had a minimum at 6 clusters ( Fig 4J). Interestingly, when whole activation patterns (spike timing over 59 electrodes) were reduced to 2 principal components, they were clustered into only two motifs, whereas the patterns reduced to 3 PCs were clustered into 6 motifs. We illustrated all patterns with a colour scale on the 2 PC plots for K-means clustering and EM 2PC clustering (Fig 4A and 4D) and on the 3 PC plot for EM 3PC clustering (Fig 4H) to visualize the clustering results. Even without clustering, the patterns in the 3 PC space were visually identified as the 6 cluster set, whereas K-means and EM 2PCs could only identify 2 clusters. However, we observed the 2 motifs identified using the other algorithms among the 6 motifs (Fig 4B, 4E and 4I). On average, the optimal number of motifs (minimum DB index) calculated using EM clustering for 2 PCs (Fig 4G) and EM for 3 PCs (Fig 4K) displayed a clear minimum value of 2 clusters (n = 6). Interestingly, in all cases, a visual inspection of the clustering identified two major motifs associated with global spike propagation pathways across the MEA. We applied the analysis described below to emphasize this finding.
The activation patterns were characterized by one major direction of a spike propagation pathway. For each pattern, we estimated the angle of the major direction by averaging all 59 vectors ( Fig 5C). Next, the clustering of the major direction angles revealed two major direction motifs in the raster plot (Fig 5A and 5B). The value of the DB index was equal to 0.08 (Fig 5D), which represents two robustly separable clusters, as shown in the histogram of the major direction angles (Fig 5F). Means of the angles from two motifs were equal to 29˚and 302˚, and the difference was statistically significant (t-test, p<0.01). Interestingly, clustering of the activation patterns composed of spike timings (EM 2 PCs, Fig 4E) and major directions angles showed similar dynamic patterns (Fig 5G and video representation S1 Video). The average DB index of the major direction angle clustering showed that the patterns clearly clustered into two major directions in all cultures (n = 6 cultures) with long superbursts (Fig 5E). The minimum value of the DB index in the cluster estimation was equal to 0.33±0.27 (mean and s.d.). Notably, a DB index with a value less than 1 indicated well-separated clusters when the inter-cluster distance was greater than the intra-cluster volume.
We shifted from analysing long superbursts to analysing regular superbursts for the presence of stable spike propagation pathways to determine whether this form of well-defined bidirectional dynamics is a unique feature of long superbursts. Regular superbursts are also composed of an initial burst (100-150 ms) followed by 3-10 small bursts (Fig 5H). The profile of the small bursts within regular superbursts was less clearly organized than that in the long superburst. In many cases, we were not able to clearly separate the small bursts due to the high variability of activity during development. In this example, 2 clusters could not be estimated correctly because one cluster contained more than 95% of all patterns. Additionally, the DB index was not monotonic (Fig 5I), in contrast to long superburst clustering (Fig 5D), suggesting the absence of the motifs. Indeed, the histogram of all major directions from this raster plot (Fig 5K) indicated the existence of one cluster (e.g., motif) with an average angle of 192˚. The characteristics of the DB index averaged over 5 raster plots (n = 5 cultures) did not indicate any clustered structure. Notably, for 2 clusters, the average DB index was equal to 1.06 ±0.65 (mean and s.d., excluding one sample with <5% in one motif), which was significantly different from the mean DB index for the long superbursts of 0.33±0.27 (n = 6 cultures, t-test, p<0.05). On average, the cluster analysis of the short bursts did not reveal a clear minimum for the DB index in 2 clusters (Fig 5J, boxplot, red lines-median values), in contrast to the recordings containing long superbursts (n = 6 cultures) (Fig 5E). Thus, the patterns in the long superbursts were clustered into two motifs associated with significantly different spike propagation directions, whereas regular superbursts did not show this feature. The mean DB index for the long superbursts 0.33±0.27 was less than 1. These low values indicated the presence of two clusters without inter-cluster overlapping [52] and, hence, are treated as statistical evidence of error-free clustering.
Then, we analysed the reproducibility of motif appearance in a long superburst sequence. Motifs in the long superburst were represented using raster plots (Fig 6A). The vertical black lines in the plot represent the motifs in the sequence of small bursts in all long superbursts. We applied EM clustering to the long superburst based on the motif frequency. This analysis identified two clusters associated with the two types of long superbursts, which are marked in blue and pink in the motif raster plot (Fig 6A, top). Each type appeared randomly in the sequence. Surprisingly, the probability of switching between the two types of long superbursts was 50% ± 13% (n = 6 cultures), indicating the random nature of this activity on a macroscopic timescale. Each long superburst comprised 100-150 small bursts, which were clearly associated with a certain major angle of the motif (Fig 6B). Surprisingly, the first type in the presented example was mostly composed of small bursts of motif #1 (96.9% of all small bursts) and, to a lesser extent, small bursts of motif #2 (3.1%). This finding revealed the presence of stable and directed spatio-temporal patterns of spike propagation pathways during the long superburst activity. In contrast, the other type of long superburst was mostly associated with motif #2 (motif #1-18%, motif #2-82%), representing a different direction of spike propagation. shows the dynamic patterns of the two estimated motifs,in which these major directions are clearly visible. On average, the probability of the appearance of each motif within its own type Blue and pink bars mark the two types of long superbursts observed in one recording. Black vertical lines illustrate bursts of a particular motif type inside each superburst. Note that one superburst type (blue) was associated with motif #1, and the other (pink) was associated with motif #2. (B) Histogram of major directions for small bursts from cultures displaying long superbursts. The two coloured clusters in the histogram represent the two motifs after clusterization, which appeared in (A). Motif dynamic patterns #1 and #2 are illustrated in C with respect to B. In all cases, the motifs from different types of superbursts have clearly different activation gradients (arrow directions) and major directions (B). The other culture showed similar principal results using the same analysis: two motifs were associated with two types of long superbursts (D), and the bursts from the motifs had significantly different major directions (E, DB index <0.1) and spatio-temporal activation patterns (F). (G) Switch probability of burst patterns during a long superburst (motif) and between subsequent long superbursts (superburst type). of the long superburst was equal to 91.5% ± 4.7% (n = 6 cultures). This remarkably high appearance of the motif in the long superbursts clearly indicated the presence of a stable functional structure of the network and, hence, reproducible dynamics.
We also measured the probability of switching between the motif types in the whole raster plot without considering long superburst indexing to quantify the stability of the motifs in the small burst sequences (Fig 6A). The switch probability was quite low, 9% ± 5% (n = 6), indicating that motifs switched quite rarely. Thus, the spontaneous bursting activity employed two basic spike propagation pathways (types of motifs) that were activated and sustained during long superbursts (10-20 sec). Representative images of similar bidirectional activity of the bursts from another culture are shown in Fig 5D, 5E and 5F. Notably, motifs from different long superburst types eventually show significantly different DB index directions of spike pattern propagation using major direction measures or dynamic pattern representations. For each motif, we estimated the mean major angle and the difference between two mean angles. On average, the difference between these angles was 95˚± 31.1˚(mean± s.d. n = 6 cultures). Surprisingly, these almost perpendicular spike propagation pathways were spontaneously selfreplicated in 6 cultures.

Discussion
In the present study, the neuronal networks formed by mature hippocampal cultures (30 DIV and older) generated specific network activity with surprisingly long sequences of bursts, i.e., long superbursts of up to hundreds of constituent bursts, with highly regular spiking patterns. In previous studies, superburst activity was reported to display a much shorter duration (up to ten bursts) [26]. In our experiments, we also observed similar activity, but in addition, more than 70% of the cultures (8 of 11) at DIV 30-35 began to spontaneously generate long superbursts with durations of up to hundreds of seconds. Six of 11 cultures generated more than 6 long superbursts during at least 30 minutes. The other cultures generated less than 2 long superbursts with a regular superburst in the background. The precise biophysical mechanism underlying the long superbursts is still largely unknown. We hypothesize that mature cultures (DIV from 35) spontaneously organize into networks with an optimal excitatory-inhibitory balance in which cycling dynamics represent a homeostatic ("natural mode") pattern that "sustains" the functional connectivity.
This long superbursting activity persisted for several days. Each superburst consisted of an initial burst with the highest firing rate, followed by a subsequence of small bursts with a 50-100 ms duration and a relatively stable 100-200 ms interburst peak interval (IBPI) (Fig  1). This interval corresponded to the bursting frequency 11.2 ± 1.5 Hz (n = 5) (Fig 2), which represented hippocampal rhythmic activity [49]. The observation of this type of periodic bioelectrical activity in a form of theta oscillations in cultures is considered a fundamental feature of hippocampal network formation, which has been widely investigated in vivo and in slices in vitro [7][8][9][10][11]. In hippocampal neural networks, the theta frequency ranges from 4-10 Hz, and the beta ranges from 10-30 Hz [53,54]. Other authors define the theta frequency as ranging from 4-12 Hz and beta frequency as ranging from 12-30 Hz [55]. Our analysis indicated a fundamental feature of hippocampal networks that generate theta rhythm, which is involved in behaviour (active motor behaviour [49]), memory tasks and the complex electrophysiological signatures of the theta and beta frequencies during the sleep state [14]. Low-frequency synchronized firing in cortical cultures is associated with classical sleep signatures [14]. The lack of external stimuli during neural network development may trigger stable low-frequency spiking activity in cortical and hippocampal cultures and indicate a functional state of sleep [14].
The spiking rate in the small bursts was 139.0 ± 78.6 Hz (n = 6), which may be associated with high-frequency oscillations and sharp wave-associated ripples in the hippocampus (100-250 Hz) [7,10] or fast gamma oscillations (90-140 Hz) [56]. A high spiking frequency has also been observed in regular superbursts, and, hence, is not directly associated with the unique characteristics of the observed long superbursts. However, the interplay between the dynamics of small bursts and the generation of high-frequency spiking plays important functional roles in vivo. High-frequency oscillations have been found to be modulated by slow theta activity in the isolated rat hippocampus [57] and in vivo [58]. Interestingly, coupling of the theta and high-frequency oscillations has been observed during rapid eye movement (REM) sleep, slowwave sleep and immobility behaviour, whereas ripples are associated with memory consolidation [10].
We postulate that this unique and rare activity appeared in cultures in which a specific balance of morphology and cell density developed spontaneously from the specific initial conditions due to mechanisms of self-organization. The cell density of neurons that survive to DIV 30 dramatically decreased from 2500 cells/mm 2 to 150 cells/mm 2 [36]. In another study, the cell density decreased from 5000 cells/mm 2 to 2500 cells/mm 2 after two weeks of development in culture [59]. Hippocampal cultures (E18) with a high density of 10 6 cells/mm 2 grown in DMEM generate superbursts with a long duration ranging from 46-91 seconds and a relatively low interburst frequency of 0.4-4.6 Hz [43], which is associated with epileptiform activity [60]. Superbursts have also been observed in E17-E18 rat hippocampal cultures with a high density of 7800 cells/mm 2 [42,61], but a detailed analysis of activity was not presented. Superbursts with a frequency of 10 Hz were only observed in the presence of increased cAMP concentrations in cultures with a density of 600 cells/mm 2 [62] or with the addition of dissociated striatal inhibitory cells [47] to cultures with a density of 1000 cells/mm 2 . In our cultures, the initial cellular density was 15000-20000 cells per mm 2 . The plated culture contained approximately 250000 cells and formed 4-5 layers. To our knowledge, our study is the first to record spiking activity under these plating conditions and stage of culture development. Thus, we propose that the plating density of the cultures was the major factor responsible for the appearance of this type of activity.
The culture medium was changed every 2 days in our experiments, whereas in other studies, the medium was changed once or twice a week [15,26]. Therefore,the frequent changes of small amounts of medium minimized the physiological stress and sustained the homeostatic balance during culture development, which may have affected the stable activity of the cultured network. The densities of glutamatergic and GABAergic synaptic terminals increased during the first 3 weeks in vitro and then saturated by . In cortical [63] and hippocampal cultures [64][65][66], the ratio of glutamatergic to GABAergic receptors in synapses and somata during development showed a similar trend to the ratio observed in vivo. The ratio of inhibitory and excitatory cells and the cell density present at DIV 30 in our experiments were also important components for the development of the in vitro neural networks with such remarkable features of activity and will be investigated in detail in further studies.
Note that we used the culture medium with a low serum concentration (only 0,4%) which is much less than usually used for glial cultures (10%) [67][68][69]. A neuron/glia ratio decreased after two weeks of culture of hippocampal cells with high density in the conditions of low serum concentration, indicating proliferation of glial cells [39]. The neuron/glia ratio did not change between 2nd and 3rd weeks of the culture [39].The inhibition of astrocyte proliferation can be explained by the effect of accumulated extracellular matrix [51,[70][71] and the contact inhibition of proliferation [72].
Theta rhythmic activity in the hippocampus has been reported to be induced and modulated by external signals originating from the entorhinal cortex. Based on our results, the spiking activity in the range of the hippocampal theta rhythm can be generated in isolated networks of hippocampal cells in the absence of stimulation under stable homeostatic conditions. Further studies using immunohistochemistry will reveal key features of cultures grown under theseconditions.
Importantly, the superbursting activity observed in cortical cultures did not exhibit the long and stable rhythmic activity reported here [14,22,26]. The reverberation of activity in the form of periodic synchronized bursts on a time scale of hundreds of milliseconds emerged or was modulated only in response to the suppression of inhibitory synaptic transmission [18,48]. Superbursting activity in dense cortical cultures (1000-5000 cells/mm 2 ) from rats and mice consisted of a sequence of bursts ranging from 0.25 to 1.25 Hz [26,46,[73][74][75]. This activity was associated with epileptic seizures in vitro [60]. In rat hippocampal slices, similar epileptiform bursts had a high amplitude (1 mV) and a low-repetition frequency (0.5-1.5 Hz), whereas theta oscillations had a low amplitude (0.5 mV) and a high frequency (5-14 Hz) [76,77]. Moreover, in postnatal cortical neurons cultured at a high density (4000 cells/mm 2 ), superbursts were generated with a frequency of 10 Hz in some cases and were induced by inhibitors of GABAergic synaptic transmission [18]. Therefore, the observation of bursting activity in the theta and delta ranges may be unique to the hippocampal cultures.
Notably, the recording area of the MEA (1.6 x 1.6 mm) was located in the centre of the circular culture and had an approximate diameter of 4-5 mm. These spatio-temporal patterns may be part of a global cycling activity with highly stable specific features of the hippocampus in vivo-theta rhythmic oscillations. The cycling pattern of the bursting activity might be triggered by pacemaker neurons [78] or may be self-organized in spiral wave dynamics. Similar spike propagation patterns have been observed in cortical cultures with chemically mediated inhibition in the network [13]. Disinhibition of GABAa-mediated synaptic transmission by bicuculline induced episodes of seizures composed of stable, short burst subsequences with 2-3 sec interburst intervals. Further studies using high-density MEA systems or fast CCD cameras for calcium imaging [60] to observe the activity of the whole culture will address these issues.
By analysing the profile of the spiking patterns in the long superbursts, we found that these patterns also become well organized and contain a small number (2-4) of basic motifs, in contrast to the regular superburst activity (Figs 4 and 5). Based on the results obtained from different clustering methods,the patterns of first spike timings in the bursts were segregated into two clusters (motifs).
These motifs defined the presence of two basic types of spike propagation direction in the burst activation pattern. These two "functional" directions further defined the activity in the form of wave-like bidirectional firing patterns that were repeated from burst to burst (Support S1 Video). The stability of the motif appearance within single long superbursts was 91.5% ± 4.7% (n = 6 cultures). This remarkable re-entry of a stable pattern clearly reveals the stable functional structure of the network. Considering the spike timing variability in the culture and clusterization inaccuracy,the motif uniqueness is likely even closer to 100%. Notably, the angle between two major spike propagation pathways of the small bursts' activation patterns was 95± 31.1˚(mean± s.d. n = 6 cultures) (Fig 6B and 6E). These almost perpendicular activity propagation pathways were spontaneously self-organized in 6 cultures and were remarkably stable during rhythmic bursting. The long superbursts were also clearly clustered into two types according to the motif appearance in the small burst subsequence (Fig 6A). Each motif appeared mostly within its own type of the long superburst, with a high probability of 91.5% ± 4.7% (n = 6 cultures), indicating that only one of two motifs was generated during each long superburst. We postulate that the dynamics of the neural network were quite stable and reproducible on timescales of milliseconds (activation patterns), seconds (small bursts) and tens of seconds (long superbursts). Interestingly, on a timescale of minutes,in which several long superburst were observed in the activity pattern,the motifs switched from one motif to another with a probability of 50% ± 13% (n = 6 cultures). The switch mainly appeared during the first small burst in the sequence. The first initiation burst with a longer duration and a higher firing rate in the long superburst sequence determined that motif. Therefore, the spiking patterns with unique orientations of activity propagation exhibited an interplay between two complex dynamic states in the network with a stochastic switch on a timescale of several minutes that initiated rhythmic activity with a precise activity pattern on a timescale of seconds and milliseconds. This conclusion complements the results of a regular superburst study [22] and can be further extended to the modelling of brain dynamics during development. Based on our results, hippocampal neuronal cultures display activity with features similar to in vivo conditions. Further studies of plating protocols and neuroengineering methods mimicking realistic hippocampal tissue conditions may identify key factors involved in the development of a functional structure in neural networks.
Notably, the well-organized global dynamics of spontaneously developing culture networks are encoded by a binary sequence and are actually represented as a telegraphic signal conveying information about the functional state of the system. We sincerely believe that this stability and reproducibility of the network states will further permit the control of the switching between the states in mature cultures and that these cultures will be useful in the design of living networks with definite functional properties in hybrid information processing systems (neurally controlled robots, "brain-on-chip", etc.).