Stimulation triggers endogenous activity patterns in cultured cortical networks

Cultures of dissociated cortical neurons represent a powerful trade-off between more realistic experimental models and abstract modeling approaches, allowing to investigate mechanisms of synchronized activity generation. These networks spontaneously alternate periods of high activity (i.e. network bursts) with periods of quiescence in a dynamic state which recalls the fluctuation of in vivo UP and DOWN states. Network bursts can also be elicited by external stimulation and their spatial propagation patterns tracked by means of multi-channel micro-electrode arrays. In this study, we used rat cortical cultures coupled to micro-electrode arrays to investigate the similarity between spontaneous and evoked activity patterns. We performed experiments by applying electrical stimulation to different network locations and demonstrated that the rank orders of electrodes during evoked and spontaneous events are remarkably similar independently from the stimulation source. We linked this result to the capability of stimulation to evoke firing in highly active and “leader” sites of the network, reliably and rapidly recruited within both spontaneous and evoked bursts. Our study provides the first evidence that spontaneous and evoked activity similarity is reliably observed also in dissociated cortical networks.

1 is also classified as ML (3.33% over the total). Therefore, we can conclude that no significant statistical difference between the firing activity of ML and follower channels should be due to a tendency for ML to record more than one neuron with respect to followers.

Properties of MLs: spontaneous activity
In our dataset we tested whether the set of MLs in a given culture remains stable across hours of  considering all electrodes of all recordings). In grey, we reported the best log-normal fit, whose equation is ).
We then analyzed the firing properties of the identified MLs. In Supplementary Figure 1b we reported the firing rate profile of each recording channel of a selected culture during 15-h recording (same culture of Supplementary Figure 1a). We depicted in red the MLs and in grey the followers. The thick black curve represents the network average firing rate. All MLs lie above the network average during the whole recording, but there are also many channels with higher firing rates which are not classified as global MLs.
Therefore, we can conclude that MLs, as we defined them, are not only apparent leaders of the bursting activity, but are also highly active.
We also measured other statistics of firing and bursting activity of MLs, namely normalized firing rate, ratio of spikes within bursts, normalized firing rate within bursts, and normalized burst duration (cf.
Supplementary Figure 1). We found that MLs feature higher global firing rates, and their activity is almost fully included within bursts, showing longer burst durations, but not higher rates within bursts (cf. Supplementary Figure 1c-f). Those results indicate that MLs are highly active and bursting neurons, and do not show tonic or persistent random spiking activity.
Additionally, inspired by the fact that also in vivo a few neurons actually feature higher firing rates than the majority of all others 1-3 , leading to a pseudo-lognormal distributions of logarithm of firing rate 4 , we derived the same probability density function of logarithm of spontaneous firing rate from our data (considering the whole dataset) and we obtained a very similar relationship. As already shown by the model hypothesized by Roxin and colleagues 4 , there is an excess of neurons firing at low rates and a lack of neurons firing at high rates compared with the best lognormal fit (cf. Supplementary Figure 1g). As also detailed in the main text, MLs show stronger responses to stimulation than followers in the early phase. This also confirms the fact that MLs tend to be highly active, both in spontaneous and evoked activity periods. Moreover, when stimulation is delivered from MLs, responses show shorter latencies (both early and late, cf. Results section in the main manuscript), indicating that the network is more promptly entrained by ML than by follower stimulation (cf. Supplementary Figure 3).

Spontaneous patterns constrain evoked patterns: multidimensional scaling results
We also used an alternative method to visualize and quantify the similarity between spontaneous and evoked patterns, which does not rely on previous clustering results. Similarly to what had been done in 5 we applied multidimensional scaling (MDS) 6 to the full matrix of all patterns' distances, in order to get a reduced 2D representation of NB patterns, given their multidimensional pairwise distances (cf.
Supplementary Methods). In fact, MDS is a nonlinear method used to represent high-dimensional datasets in a low-dimensional (typically 2D) space such that pairwise distances are preserved as well as possible (i.e. points which are close in the original high-dimensional space will also be placed close by in the 2D projection) 5 .

Spike sorting
The results presented in the main manuscript were obtained by applying no spike sorting procedure.
This choice was made according to many other published studies 7-10 , also reporting that in such cultured were saved for further analysis. All spikes were aligned to their minimum at data point 10. In order to avoid spike misalignments due to low sampling, spike minima were determined from interpolated waveforms of 64 samples, using cubic splines. Channels were considered for sorting if recording at least 50 spikes during the whole experiment. First, unit clusters were determined by analyzing the first 30-minute recordings (and considering no more than 20,000 spikes). All other spikes detected during the recording were assigned via template matching, i.e. they were assigned to the closest cluster unless they were too far (threshold: 3 times standard deviation from the cluster center). The obtained results were further refined and confirmed by visual inspection of an expert user. Unclustered spikes were not considered.

Normalization of pattern distance
Normalized distances following the approach described in the main manuscript (and illustrated in Supplementary Figure 5a) do not depend on pattern length and, when re-ordered according to the pattern leader, nicely cluster along the diagonal. In Supplementary Figure 5b we reported a representative example of the effect of normalization on pattern distance values. Upper panels show raw distance matrices, whereas lower panels depict normalized distances, re-ordered either according to pattern length (left) or to pattern leader (right). It is evident that raw distances strongly depend on pattern length, preventing clusters of similar patterns to be detected, whereas normalized distances do not.

Additional information on pattern clustering
According to 12 , the distance matrix is first rearranged by a standard agglomerative dendrogram method, based on Euclidean distance (provided by the Statistics Toolbox in Matlab). Second, the algorithm iteratively looks for sets of consecutive N patterns in the re-ordered matrix that satisfy two conditions: separation, i.e. a pattern belonging to a given set will have the lowest average distance with patterns belonging to the same set than to any other set, and isolation, i.e. no pattern will belong to more than one set. The parameter N is chosen equal to the square root of the total number of NBs included in the analysis, as suggested in the original paper 12 . Hence, the k selected sets (i.e. clusters) do not overlap and they maximize the inner similarity (i.e. minimize the inner distances among patterns) (cf. Supplementary Figure   6). Once the number of different clusters is identified, they are used as templates and NBs not included in any cluster are associated via template matching to the cluster they are closer to (i.e. a NB is included in a cluster if its average distance from the patterns belonging to that cluster is less than the average plus three times the standard deviation of that cluster's inner distances). The clustering procedure is applied separately to both spontaneous and evoked activity. In case of evoked patterns, the procedure is applied independently to each stimulated channel's responses, and only one cluster per channel is considered for further analysis (the one including the highest number of responses

Multi-dimensional scaling (MDS)
Non-classical metric multidimensional scaling (MDS) was performed in MATLAB 6 . For this analysis we considered up to a maximum of 2500 spontaneous patterns and all evoked patterns, as also done for previous analyses on pattern similarity. MDS was applied to the full matrix of normalized edit distances between all pairs of patterns (spontaneous and evoked). We plotted the results of MDS in two dimensions, and we used the convhull function (provided by Matlab, useful to determine the convex hull of a set of points in 2-D space) to find the convex polygons delimiting either spontaneous or evoked patterns. We then used the inpolygon function (also provided by Matlab) to determine the ratio of evoked patterns included in the hull defined by spontaneous patterns (separately for each experiment). We then transformed our plot into a binary image (using poly2mask, Image Processing Toolbox, Matlab) where pixels included within hulls are non-zero, to compute the intersection area between hulls for each experiment.