Multiscale functional connectivity estimation on low-density neuronal cultures recorded by high-density CMOS Micro Electrode Arrays

Journal of Neuroscience Methods 207 (2012) 161–171

Contents lists available at SciVerse ScienceDirect

Journal of Neuroscience Methods journal homepage: www.elsevier.com/locate/jneumeth

Computational Neuroscience

Multiscale functional connectivity estimation on low-density neuronal cultures recorded by high-density CMOS Micro Electrode Arrays Alessandro Maccione a , Matteo Garofalo a,b , Thierry Nieus a , Mariateresa Tedesco b , Luca Berdondini a , Sergio Martinoia a,b,∗ a b

Department of Neuroscience and Brain Technologies - Istituto Italiano di Tecnologia, via Morego, 30, 16163 Genova, Italy Neuroengineering and Bio-nano Technology Group (NBT), Department of Biophysical and Electronic Engineering (DIBE), University of Genova, Genova, Italy

a r t i c l e

i n f o

Article history: Received 5 December 2011 Received in revised form 30 March 2012 Accepted 2 April 2012 Keywords: Functional connectivity Effective connectivity CMOS Micro Electrode Arrays Low density cultured neurons

a b s t r a c t We used electrophysiological signals recorded by CMOS Micro Electrode Arrays (MEAs) at high spatial resolution to estimate the functional-effective connectivity of sparse hippocampal neuronal networks in vitro by applying a cross-correlation (CC) based method and ad hoc developed spatio-temporal filtering. Low-density cultures were recorded by a recently introduced CMOS-MEA device providing simultaneous multi-site acquisition at high-spatial (21 ␮m inter-electrode separation) as well as high-temporal resolution (8 kHz per channel). The method is applied to estimate functional connections in different cultures and it is refined by applying spatio-temporal filters that allow pruning of those functional connections not compatible with signal propagation. This approach permits to discriminate between possible causal influence and spurious co-activation, and to obtain detailed maps down to cellular resolution. Further, a thorough analysis of the links strength and time delays (i.e., amplitude and peak position of the CC function) allows characterizing the inferred interconnected networks and supports a possible discrimination of fast mono-synaptic propagations, and slow poly-synaptic pathways. By focusing on specific regions of interest we could observe and analyze microcircuits involving connections among a few cells. Finally, the use of the high-density MEA with low density cultures analyzed with the proposed approach enables to compare the inferred effective links with the network structure obtained by staining procedures. © 2012 Elsevier B.V. All rights reserved.

1. Introduction The importance of estimating structural, functional and effective connectivity at different scale of complexity (i.e., from a few cells to organized tissue or whole brain) has been recently addressed by many works and it is gaining momentum within the neuroscientific community. The concept of connectome is now well established (Sporns et al., 2005) and it has individuate a new field of study in which investigations about the interplay between dynamics and connectivity plays a major role. The increasing attention about functional interaction among neurons arise from the widely accepted hypothesis that major brain functions are executed through the joint actions of neurons or cell assemblies. In fact, it has been shown in many reports, and specifically in brain–machine interface studies (Nicolelis and Lebedev, 2009) that information is not encoded by single cells but rather by populations (Bi and Poo, 2001). From this perspective, it becomes fundamental to estimate functional and effective interconnections

∗ Corresponding author at: Department of Neuroscience and Brain Technologies - Istituto Italiano di Tecnologia, via Morego, 30, 16163 Genova, Italy. E-mail address: [email protected] (S. Martinoia). 0165-0270/$ – see front matter © 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.jneumeth.2012.04.002

among neurons and possibly to identify directionality and synaptic pathways. Thus, reliable methods to investigate and estimate how single or groups of neurons are connected (Garofalo et al., 2009; Rieke et al., 1997) becomes essential as dynamically interacting neuronal populations represent functional blocks on which complex behavior such as learning or memory are based (Harris, 2005). In general, statistical or information theory based methods, applied to the recorded neuronal signals, allow to localize information flow and to build maps of the relevant functional-effective connections. In a recent review by Feldt et al. (2011) the issue of structure and function is addressed by underlying that to bridge the gap between cellular and behavioral levels will require an understanding of the functional organization of the underlying neuronal circuits. One possibility to unravel the complexity of neuronal networks is therefore to understand how their connectivity emerges during development. To this end in vitro studies (as specifically shown also in the cited review) can contribute to lay the foundations to the interplay between structure and function and contribute in answering some of the central questions related to the interplay among dynamics, functional connectivity and network morphology. In this work, we present a specific application of a crosscorrelation based method with additional filtering procedures, to

162

A. Maccione et al. / Journal of Neuroscience Methods 207 (2012) 161–171

study functional-effective connectivity inferred by electrophysiological recordings in hippocampal cultures. To this end we used a recently introduced high-density MEAs (Maccione et al., 2010) (commercially available from 3Brain, Switzerland) coupled to low-density cultures to better access single cell activity recordings together with the possibility to visualize interconnected networks. These measurements would have been unfeasible with conventional MEA systems characterized by good sampling rates (i.e., up to 50 kHz), but with low number of recording sites and low spatial resolution (i.e., typically 60–256 microelectrodes of 10–30 ␮m in diameter). The CMOS-MEA chip we used is based on the Active Pixel Sensors (APS) technology (Berdondini et al., 2009a) enabling acquisitions from 4096 microelectrodes at a full frame rate of 7702 Hz for each channel and at a spatial resolution of 21 ␮m. Thus, the use of APS-MEA devices, combined to sparse cultures, allows to map neuronal signaling in large-scale networks at spatial resolution down to the cellular level up to a possible identification of its anatomical connections (i.e., the set of physical or structural-synaptic connections linking neuronal units at a given time; Sporns and Tononi, 2002). In this article we focus on the estimation of the functionaleffective connectivity from extracellular electrophysiological recordings by applying cross-correlation based methods (Salinas and Sejnowski, 2001) and additional spatio-temporal filtering procedures. The rationale behind this approach relies in the possibility to infer effective links related to possible causal relationships between pairs of neurons estimated from the statistical analysis of the spontaneous activity. The connection strength (correlated to the synaptic weight) between two neurons is considered to be proportional to the value given by the method. By using the above introduced experimental platform we analyzed four low-density cultures (see Table 1) to estimate their functional-effective links and to reconstruct connectivity maps superimposed to fluorescent morphological images of the cultures. All the analyses are based on CC methods with additional filtering procedures that are used to discriminate between possible causal relationships and spurious connections, and thus to improve the reliability of the estimated maps. These results have been achieved by combining low density cultures with high density MEAs and they illustrate the potentialities of the approach (including immunostaining) and highlight relevant future applications for the high-density APS-MEAs.

2. Materials and methods 2.1. Cell cultures and experimental set-up Hippocampal neurons were dissociated from rat embryos and plated on APS-MEAs previously coated with adhesion promoting molecules (poly-d-lysine and laminin) as already presented in the literature (Berdondini et al., 2009a). The nominal concentration varies from 1 to 4 × 103 cells/device (cf., Table 1) resulting in a final density ranging from about 80 to 200 cells/mm2 . After plating on the chip, cultures were filled with 1 ml of nutrient medium (Brewer et al., 1993) and placed in a humidified incubator (5% CO2 )

at 37 ◦ C. After about two weeks in vitro cultures start showing a spontaneous electrical activity composed by isolated spikes and short bursts involving many channels (neurons) within the whole network. All measurements started 20–30 min after positioning the dishes upon the recording device. On the time scale of our measurements (up to 45 min) the recordings were rather stable in terms of their spiking activity (data not shown). The experimental platform, based on the Active Pixel Sensor (APS) MEA, was already presented in (Imfeld et al., 2007, 2008). Briefly this innovative system enables multisite recordings from 4096 electrodes arranged in a 64 × 64 layout. Each squared electrode of 21 ␮m side is spaced at 42 ␮m (pitch) and sampled at 7702 Hz (Fig. 1D provides typical raw data traces recorded from three different electrodes). Custom software was developed to manage the large throughput of data (about 60 MB/s) allowing real-time signal visualization and processing based on image/video analysis. The recorded signals (an example in Fig. 1E), ranging from random spiking activity to more complicated and synchronized burst behavior, were detected by using the Precise Timing Spike Detection (PTSD) algorithm (Maccione et al., 2009) that allows a fast and precise temporal identification of the spikes. Bursting activity was detected by using a single channel spike train burst detector as presented in previous paper (Chiappalone et al., 2005). Additionally, first order statistical parameters such as the Inter-Spike-Interval (ISI), the Inter-Burst-Interval (IBI), the Mean Firing Rate (MFR) and the Mean Bursting Rate (MBR) (Rieke et al., 1997) were computed (results not shown). No sorting methods were performed since electrodes record signals mainly from a single neuronal unit by using low density culture plated on the high density chip as demonstrated in a previous paper (Berdondini et al., 2009b).

2.2. Staining procedures After spontaneous electrophysiological recordings, low density cultures were fixed and stained. A morphological and topological representation of the entire network was obtained by collecting and then stitching several pictures at a sufficient high resolution to identify neuronal processes. These images were then used to present functional connectivity maps by superimposing the estimated links to the images of the network obtained by the staining procedures. Briefly, as already described in previous papers (Berdondini et al., 2009a; Maccione et al., 2010), cells on the chip were fixed at room temperature at 4% paraformaldehyde in PBS, permeabilized for 10 min and exposed for 30–40 min to a blocking buffer (in some cases we were able to reuse the devices after staining). After 2 h incubation cultures were labeled with the secondary antibody Alexa-Fluor488 conjugated goat-antirabbit for MAP2. Then for culture hn3 (see Table 1) we also labeled with AlexaFluor 546 conjugated goat-antimouse for NeuN (see Fig. 1B), while for cultures hn2 and hn4 we used DAPI (Sigma–Aldrich). Afterwards cultures on the chip were inspected under BX51M Olympus microscope equipped with 20 × 0.45 SLM objective. DIC and fluorescence images were taken with the care of preserving overlapping regions

Table 1 Statistical parameters of the analyzed cultures. The first 2 columns provide the nominal concentration and the dimension of the drop used to seed cells, DIV (Days In Vitro) column refers to the age of the cultures when measured. MFR and MBR are calculated on the whole duration of the recording. Last column relates to the number of links after the thresholding procedure (cf. Section 2). Topological parameters

hn1 hn2 hn3 hn4

Statistical analysis

Cell/␮l

␮l

Imaging

DIV

Duration (s)

MFR (spike/s)

MBR (burst/min)

Links

∼190 ∼240 ∼120 ∼170

20 12 11 25

No Yes Yes Yes

28 19 33 18

2700 1200 1200 300

0.42 0.69 0.37 0.48

8.4 11.9 7.85 5.73

1928 1744 205 169

A. Maccione et al. / Journal of Neuroscience Methods 207 (2012) 161–171

163

Fig. 1. High-density CMOS device coupled to low-density hippocampal cultures. (A) Packaged chip; (B) active area with low density culture stained with MAP2 (green) and NeuN (red); (C) zoom on single pixel-single neuron. (D) Raw data traces recorded by the two electrodes under isolated neurons in (C). Synchronous activity made of isolated spikes is clearly visible. (E) Example of raster plot (about 200 active electrodes) showing the global dynamic behavior of the network. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

by using a Hamamatsu ORCA ER II camera driven by Image-Pro Plus software (from Media Cybernetics). Successively images were both manually stitched together by using Photoshop CS3 (from Adobe) or automatically stitched by using open source software tool as ImageJ (http://rsbweb.nih.gov/ij/). The result is a high-resolution picture of the entire active area of 2.7 mm × 2.7 mm, where neuron’s nuclei and neurites are clearly resolved. As an example in Fig. 1B and C nuclei and neurites growing on the electrode array are clearly visible. 2.3. Cross-correlation Cross-correlation (CC) functions (Salinas and Sejnowski, 2001) were evaluated among the networks’ nodes (i.e., electrodes) by considering each possible pair of spike trains extracted by each active electrode (i.e., firing rate > 0.1 spike/s) of the CMOS-MEA. A Cross-Correlogram is a temporal function which combines the

firing information of one target neuron with a reference one. Considering the low density of the network that allows identification of the sources, CC was used to detect possible causal interactions between neuron’s pair. The connection strength is estimated on the basis of the peak value of the CC function (Cpeak ) and the directionality is derived from the temporal position of the corresponding peak latency (i.e., before or after t = 0). By considering the peak amplitude of each CC function, a Connectivity Matrix (CM) is defined where the highest values are supposed to correspond to the strongest connections (Garofalo et al., 2009 Plos One). The method we used for computing the Cross-Correlogram Cxy () is a mixture between the classical definition (Knox, 1981; Rieke et al., 1997) and the approach proposed by Eytan and Marom (Eytan et al., 2004). We counted the number of spikes in the Y train within a time frame (±T) around the spikes of the X train, using bins (steps) of amplitude  (set at a multiple of the sampling frequency). The Cxy () is then normalized, by dividing each element

164

A. Maccione et al. / Journal of Neuroscience Methods 207 (2012) 161–171

Fig. 2. Cross-correlation and spatio temporal filters. (A) Example of cross-correlation Gaussian fitting function (red line) derived from Cross-Correlogram (black line) with a clear peak around 1 ms. (B) Strongest 100 links delay–distance relationship on culture hn3. Other than the application of the PF (dashed black line) in order to discard non physiological delay–distance relationships, STF (gray line) identifies two different area: STF 0–1 encompassing all the physiological plausible links below the shifted time–space filter at 1 ms and STF 1 for all the link above the filter. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

of the array by the geometric mean of the number of spikes in the X and the Y train, yielding the CC function:

CXY () =



1 NX NY

NX 



(+(/2))

X(ts )Y (ts + ti )

(1)

s=1 ti =(−(/2))

where ts (ti ) indicates the timing of a spike in the X (Y) train and NX (NY ) is the total number of spikes in the X (Y) train. In this paper, CC function was evaluated on a time window T of ±150 ms and by setting a time bin of 3 frames (∼0.39 ms). In order to detect relatively fast interactions among electrode pairs, we look for the peak of the function in a sub-window of ±5 ms. Direction of the links was then accounted by looking at the latency of Cpeak . Absolute CC peaks with zero delay are discarded to ignore co-activation unless they refer to possible unresolved intra-neuron propagation. These links are identified when Cpeak falls in t = 0 and the link length is shorter than 84 ␮m; that is 2 electrodes separated at most by one other electrode and thus a distance compatible with neuritic propagation (i.e., assuming a propagation speed of 400 mm/s, see also Section 2.5). Moreover, in order to take into account only relevant CC, we varied empirically the constraint on the minimum number of events contributing to the Cpeak . We looked at a threshold able to detect robustly the links, by avoiding the presence of a high number of false positive links in the connectivity map. The resulting value, 7 occurrences for Cpeak , was set as the default parameter used to obtain all the results reported in the following. Importantly, according to these criteria almost no CC peaks could be detected outside the ±5 ms time interval, further motivating the choice of limiting the search of the peak in a restricted time interval. Thus we focused our attention to relatively fast interactions and possible causal relationships trying to improve the reliability of the estimated maps. To this end, a qualitative comparison with fluorescent images (i.e., maps superimposed to the images) is reported. In Fig. 2A an illustrative example of a significant cross-correlation function is shown. To efficiently manage the large amount of data, the CC function was coded in C# and a typical run (5 min data) lasted 2 min long.

2.4. Thresholding procedure Since the CC function (Eq. (1)) is computed across all electrodes of the MEA recording system we need to establish a criterion to define reliable connections based on their Cpeak values. The statistical significance of the estimated connections, i.e., above a certain noise level, was assessed by a shuffling procedure. We shuffled the ISIs of the spike trains maintaining their original (ISI) distribution as we aimed at destroying any causal effects, estimating the noisy Cpeak value expected between two random spike trains having the same ISI distribution and firing rate of the original ones. The shuffling procedure was applied to the spike trains yielding the highest Cpeak values and for each of them the corresponding spike trains (Xi ,Yi ) were shuffled 5000 times. Then, a shuffled CC threshold (THsh ) was defined as m + 3std (Grun and Rotter, 2010) where m is the mean value and std the standard deviation of the shuffled dataset. To be more conservative and in order to highlight the most relevant links, we then selected the strongest K-links among the possible N-links above the defined THsh . This was done in agreement with a previous presented work (Garofalo et al., 2009) in which we considered the connection strength proportional to the Cpeak amplitude and we selected the strongest excitatory links. It should be underlined that other selection criteria could be used (e.g., on the basis of percentage or percentile values) as any further selection procedure is arbitrary and it aims only at further discriminating the most relevant Cpeak values among the statistically significant ones. 2.5. Spatio-temporal filtering procedure We introduced specific spatio-temporal filters (STF) by defining a distance-dependent latency threshold (DdLT) in order to refine the significant links previously selected based only on the links strength (Cpeak ). To this aim, we computed the link length (evaluated as the Euclidean distance) and the Cpeak delay for each electrodes pair. The DdLT is defined by assuming a maximum propagation velocity (i.e., 400 mm/s), a value higher than those reported in literature (i.e., 350 mm/s) by using similar biological preparations (Bonifazi et al., 2005). If the Cpeak delay is smaller than the one allowed by the maximum propagation velocity, the link

A. Maccione et al. / Journal of Neuroscience Methods 207 (2012) 161–171

165

Fig. 3. Reliability of the detected links. (A) Stability of the functional links estimated by means of the CC respect to the duration of the observation window (chunk duration). The overlap percentage (y-axis) represents the number of connections commonly identified on different chunks of the same recording session (hn1 – 45 min). Chunks of 3, 7 and 15 min were considered. (B–D) Histogram of the CC values (Cpeak ): the threshold evaluated by means of the shuffling procedure is represented by the dashed line. Gray bars show the histogram of all the Cpeak while the black ones show the strongest 200 links for (B) hn2 and 100 links for (C) hn3 and (D) hn4.

is discarded as it is considered not physiologically plausible. We called this distance-dependent latency threshold: PhysiologicalFilter (PF, black dashed line in Fig. 2B). As an example, black crosses in Fig. 2B represent the links which do not satisfy the DdLT conditions imposed by the PF and, therefore, are discarded. Additionally, we introduced shifted spatio temporal linear filters to account for the delayed generation of spikes in the postsynaptic cell. As a reference for selecting the faster connections we adopted the fast excitatory AMPA rise time response (∼1 ms) to define the minimum acceptable delay of our filters. Thanks to their application, it is possible to select those links characterized by a certain distance-dependent latency relationship. STFs application allows us to focus possibly on mono-synaptic or poly-synaptic connections. In the following, we refer to these filters indicating the time shift introduced. As an example, STF 1 means that we are focusing on those links identified by a filter with the same slope of the PF but shifted at 1 ms. In this case, STF 1 (e.g., gray line in Fig. 2B) identifies a subset (i.e., black scattered points) of those links selected by the PF (i.e., black and gray scattered points). Finally, it is indeed possible to select refined distance-dependent latency connections by combining 2 STF filters (cf. Fig. 2B for STF 0–1). 3. Results Spontaneous activity of four low-density hippocampal cultures (hn1, hn2, hn3, hn4 in Table 1) was recorded at different Days In Vitro (DIVs) and consequently analyzed. Table 1 provides both parameters that influence network topology as cell plating density or DIVs of the experiment and main statistical features such as Mean Firing Rate (MFR) and Mean Bursting Rate (MBR) or the

number of significant estimated functional-effective links for each culture. Particularly, the longest recording was performed on hn1 and this dataset was used to evaluate the temporal stability of the estimated connectivity (see Fig. 3A). All results were obtained by determining significant functional connections (FCs) based on the following criteria: (i) FCs were computed across ‘active’ electrodes, i.e., satisfying the requirement that MFR >0.1 Hz (Section 2.3); (ii) Cpeaks had to cross two thresholds, namely (ii-a) the 7 occurrences in the non normalized CC function (Fig. 2A, Section 2.3) and (ii-b) the shuffled THsh value in the normalized CC function (Section 2.4); (iii) the FCs had to pass a Physiological Filter (Fig. 2B, Section 2.5) Finally, the obtained FCs were further restricted to the first 100–200 connections, a number that tends to minimize the links discarded by the PF (cf. Fig. 4A). 3.1. Analysis of the connectivity maps: assessment of the reliability of the CC method In order to test the reliability of the detected links as a function of the recording time duration, we analyzed how the resulting connectivity maps are influenced by changes in the temporal length of the recorded data. At first the stability of the longest available recording (hn1 – see Table 1) was assessed by verifying that the MFR did not show any particular trend in time (i.e., the slope of the MFR versus time relationship was statistically compatible with the null value). Then we divided the 45 min into chunks of 3 min and, for each segment, we computed the corresponding Connectivity Matrix. The resulting connectivity maps were compared to evaluate the percentage of connections commonly identified in all the recording chunks. This analysis was repeated by varying both the recording chunks duration (7 and 15 min) and the number of

166

A. Maccione et al. / Journal of Neuroscience Methods 207 (2012) 161–171

Fig. 4. Analysis of the connectivity parameters. (A) Physiological-filter (PF) effect on hn2: percentage of the points discarded by the spatio-temporal filter. (B) Links length of all links (gray bars) and of the strongest 200 ones (black bars). (C) Distribution of the Cpeak delays. (D) Relationship between delay and length of the strongest 200 links: black crosses represent the 17 connections discarded by the physiological-filter (black line).

connections (cf. Section 2.4). Fig. 3A shows the percentage of links commonly identified across the chunks by the CC (y-axis) as a function of the number of considered links (x-axis). The region of maximum overlap (between 30 and 50 connections in Fig. 3A) gives an indication about the number of detected links (via the CC method) that are temporally stable. Note also that, in this range, 15 min lasting chunks provide an overlap percentage always over 60% with a peak of about 80% and that 7 min chunks still perform rather well. As a whole, the presented data show that the range between 20 and 100 links always gives a reasonable overlapping among the detected links regardless of the particular recording length. The additional choice we did on the statistical relevant links detected by the implemented CC method (i.e., links with Cpeak > THsh ) was to account only for a restricted number of strongest connections. Therefore, regardless of the number of estimated links by the shuffling procedure presented in Section 2 (cf. Table 1 last columns), to be more conservative, we chose the strongest 100 connections for hn1, hn3 and hn4 and 200 ones for hn2. The choice of 200 for culture hn2 is due to the presence of a fairly big cluster of cells that concentrates much of the strongest links in that particular region (see close-up image of the cluster referring to hn2 in Fig. 6B). Following this restrictive choice, Fig. 3B–D shows the Cpeak distribution histogram of the strongest 200 (100) selected links – black bars – and the THsh value – dashed line – for hn2 (hn3 and hn4).

3.2. Functional maps and spatio-temporal selection Experiment hn2 is used to illustrate the effect of the spatiotemporal filtering procedures (i.e., PF and STF in Fig. 2B) and to characterize the estimated functionally connections. To this end, the distribution of the Cpeak delays, the histogram of the links length and their relationships were investigated. In the previous paragraph we showed the Cpeak histogram distribution (Fig. 3B–D) for three different experiments and how the adopted criterion for the threshold selection (i.e., a given number of the strongest links) suitably satisfies the lower bound condition imposed by the shuffling procedure (e.g., dotted vertical lines in Fig. 3B–D). Fig. 4A shows the percentage of connections discarded by applying the PF as a function of the considered links (up to the maximum number of significant links = 1744). With the exception of the first connections (1/40 links), a plateau at about 9% can be observed up to 300 links. Afterwards, the percentage of removed links increases steeply, demonstrating that in general many spurious connections could be identified without the PF application. Note that all these spurious connections are statistically significant on the basis of the THsh value but, just restricting the criterion to the first 200 strongest connections, many of them are discarded before the application of PF. Fig. 4D reports the space–time relationship of the 200 strongest links of the experiment hn2. Those represented by black crosses below the PF (i.e., black line) are filtered out while gray dots (i.e.,

A. Maccione et al. / Journal of Neuroscience Methods 207 (2012) 161–171

points above the line) represent physiologically plausible connections. We deeper characterized the inferred maps of hn2 by analyzing the distribution of the link length (Fig. 4B), delays (4C) and the relationship between delay and connections length with respect to the application of the PF filter (4D). Black bars refer to the strongest 200 links while gray bars to all the connections above the THsh (i.e., 1744 links). Fig. 4B highlights how the map of the strongest 200 values (black bars) is characterized by many short links and a few long ones, as it could be expected from 2D cultured networks (Shefi et al., 2002; Takahashi et al., 2010). On the other hand, by considering all the links (i.e., gray bars), the length distribution becomes flatter as we introduced many weak and indirect links into the connectivity map. Fig. 4C reports the distribution of the Cpeak delays for the first few ms as we restricted our study to the faster and possibly causal interactions. If selecting the strongest 200 connections (i.e., black bars) the links characterized by the shortest delay (equal to 1 bin; 0.39 ms) are the most detected while no preference are observed among the other delays. The same analysis, performed on hn3 and hn4 (data not shown), confirmed that a peak occurs in correspondence of the shortest delay (i.e., 1 bin) and no other significant trends were observed. On the other hand, when considering all

167

the statistically significant links (i.e., gray bars; Fig. 4C), no specific trend is observed. Finally, the scatter plot of delays vs link lengths (see Fig. 4D) shows that the faster detected connections are broadly distributed and only the application of the PF filter discard the unreasonable fast interactions for large distances. We also tested the relationship between the Cpeak amplitude and the Cpeak delay (data not shown) and we did not find any clear trend, suggesting that there is no specific dependency of the supposed synaptic strength (Cpeak ) on the time to which synaptic mediated connections can occur.

3.3. Estimated functional-effective links In this section we report the functional-effective connectivity maps of hn2, hn3 and hn4 (cf. Figs. 5 and 6) superimposed to the immuno-fluorescent image of the neuronal network. The links are represented by arrows that give indication of the direction, and are color coded proportionally to the connection strength (i.e., from white – weak to red – strong). Further, white squares identify the electrodes which are significantly cross-correlated (i.e., the Cpeak > THsh ) with at least one electrode.

Fig. 5. Connectivity maps superimposed to the corresponding morphological fluorescent image of the culture hn3. A color-code (from white – weak to red – strong) shows the links strength. (A) CC map of the strongest 100 links. (B–D) Effect of the filtering procedure. (B) Physiological-filtering (PF): links which do not satisfy the minimum delay–distance condition are considered as not physiological and discarded. (C) STF 0–1 aimed to individuate fast propagation (monosynaptic). (D) STF 1 aimed to identify poly-synaptic propagations.

168

A. Maccione et al. / Journal of Neuroscience Methods 207 (2012) 161–171

Fig. 6. Connectivity map of hn2 and hn4. 200 (100) strongest links are filtered and superimposed to the fluorescent image of hn2 (hn4) (A) Physiological-filtered map of hn2 – 183 connections; (B) closed up on the cluster of cells containing 118 of the strongest connections; (C) STF 0–1 filtered map of hn4 – 31 connections; (D) STF 1–3 filtered map of hn4 – 46 connections.

Fig. 5 refers to experiment hn3. We applied the PF (dashed black line in Fig. 2B) and the STF1 (gray line in Fig. 2B) on the strongest 100 connections. Their combination allows to individuate three regions: the PF area (i.e., scattered points – gray and black), the STF1 area (i.e., black points) and the STF 0–1 area (i.e., gray points). The links discarded by the PF, because not compatible with signal propagation speed are indicated by the black crosses. Maps relative to these three regions are reported in the illustrative example of Fig. 5. Fig. 5A shows the map estimated by selecting the strongest 100 links without using any distance-dependent latency threshold (DdLT). Fig. 5B reports the map obtained by selecting the connections (within the ones presented in Fig. 5A) belonging to the specific spatio-temporal region identified by the PF (i.e., black and gray points of Fig. 2B). The map (Fig. 5B) obtained by PF application, shows that some of the links identified by the CC disappeared (20 links are discarded, i.e., the black crosses in Fig. 2B) and as expected, among the discarded, the longest ones represent the majority. The identified 80 links are then divided into two complementary regions by the spatio-temporal filtering procedure (cf. Section 2.5) set at 1 ms (i.e., gray line in Fig. 2B). Particularly, Fig. 5C shows the 53 links identified by the STF 0–1 (i.e., gray points of Fig. 2B). These connections might be considered as the ones generated by signals propagating very fast (and efficiently) between two neurons (e.g., mono-synaptic fast connection). The shortness of the links is

a good indication, together with visual comparison with the morphological image, that possibly direct synaptic connections were identified (see also magnifications reported in Fig. 7B). Further analysis related to neuron’s pair interaction is presented in a separate section. Fig. 5D shows the links already identified by the STF1 and corresponding to the 26 black points of Fig. 2B. It should be underlined that these links could be interpreted as mono and poly-synaptic connections in which neuron pairs or microcircuits are formed to form specific assemblies. Similar as hn3 results relative to hn2 and hn4 are reported in Fig. 6. The peculiarity of hn2 experiment consists in the big cluster of neurons clearly visible in the bottom part of the image (Fig. 6A). The cluster shows a strong and synchronized electrical activity (data not shown) that is captured by the CC analysis which predicts almost all the possible connections in this area. By applying, the PF to the strongest 200 links, the map of Fig. 6A shows 183 links, resulting in a loss of almost 9% of the initially identified links (as indicated also by Fig. 4A). Of the 183 links, 118 are confined in the big cluster. It should be noted that the method is able to identify cells interconnections in the cluster but it is not able to discriminate the directionality. In fact, the temporal resolution (i.e., sampling interval of 0.13 ms and bin size of ∼0.39 ms) cannot discriminate positive or negative peak delays: the Cpeak always

A. Maccione et al. / Journal of Neuroscience Methods 207 (2012) 161–171

169

Fig. 7. Analysis of structures constituted by two and three interconnected neurons. (A) Magnification of hn2 and (B) hn3 CC maps to individuate small group of neurons after the application of the PF (C) CC are divided into 4 pools on the basis of the position of the Cpeak . Each CC pool is fitted with a Gaussian curve and its variance is evaluated. (D) Cpeak delays of the “triangle structures” of hn3 (black cross), hn4 (gray scattered points) and hn2 (black scattered points). The direct link delay (x-axis) is usually smaller than the indirect path one (y-axis).

falls in the central bin accounting for strong correlation with no information on possible causal interactions. In Fig. 6C and D, as a further illustrative example, the maps of hn4 are presented. The immuno-staining shows a very low number of neurons and connections (see Table 1). For this reason we selected only the strongest 100 links, 93 of them were maintained after the PF procedure. Maps reported in Fig. 6C and D were obtained by applying the STF at 1 and 3 ms identifying 31 (STF 0–1) and 46 (STF 1–3) connections. In this particular experiment the confinement of sub-networks in different separated areas (as it can be observed from the image) well matches the estimated connectivity maps. As shown in the illustrative examples presented above, the possible comparison between effective and structural–topological connectivity of the culture could provide, if properly elaborated, an interesting description of functional interactions among neurons at the whole network level in terms of propagating pathways and main synaptic connectivity.

3.4. Neuron pairs and microcircuit analysis Scaling down from the whole network to few cellular units we identified and analyzed microcircuits and neuron pairs. As an illustrative example, Fig. 7A and B report two details of the connectivity maps already shown in Figs. 6A and 5A respectively.

From a mathematical point of view CC detects direct as well as indirect connections since it is not able to recognize the direct causal interactions as other methods do (e.g., partial CC, Eichler et al., 2003; redundancy, Bettencourt et al., 2007). Therefore, maps inferred by CC are affected by possible multiple indirect connections. In order to investigate if indirect links can provide “false positive” connections, we focused on structures constituted by two and three interconnected neurons. Indeed, the PF filtering procedure discards some of the false positive and in the analysis of the links we started from the connectivity maps after the application of the PF filter. The analysis of the connections between neuron pairs was motivated by the observation that the bell-shaped CC had a lower spread when the time lag fell close to zero, providing an apparent relationship between the spread of the curve and its peak position. Indeed the fitting of the CC distributions with a Gaussian function shows a clear relationship between the STF delay and the CC variance (Fig. 7C). Notably, the CC curves in zero has a variance comparable to that of the autocorrelation ones, meaning that they are compatible with inter-neuronal signal propagations, while for non-zero delays the CC variance increases, showing that those estimated links are sustained by synaptic propagation. Error bars were estimated by means of a bootstrap procedure (Duda et al., 2000). The higher variances at increasing STF values can be explained by considering the higher variability of propagating signals in slow

170

A. Maccione et al. / Journal of Neuroscience Methods 207 (2012) 161–171

mono or poly-synaptic pathways compared to fast single synaptic connections. Finally, Fig. 7D reports the analysis on microcircuits involving three nodes each of them connected to the other two. For each microcircuit the delays of both the direct (delay13 ) and the indirect path (delay12 + delay23 ) were calculated. The delay of the indirect link is, in the 85% for all the three experiments here considered, bigger than the delay of the direct one. This result suggests that the direct link seems not to be dependent from the indirect path (i.e., through node 2), showing a causal connection that is independent of the indirect pathway.

4. Discussion and conclusions The issue of investigating the relation between structural and functional brain networks is one of the central topics of modern neuroscience. The problem has been widely addressed at different scales (Bullmore and Sporns, 2009) and it has been shown that there are some general topological features that are shared (and preserved) at different complexity levels ranging from whole brain scale (in human) down to cellular scale (in animals). Small word connectivity seems to be a common signature observed (as extremes) in in vitro cultures and in brain regions (Pajevic and Plenz, 2009; Sporns, 2006). Within this context, graph theory and statistical methods have been widely applied to investigate the relationships between structure and function. As recently underlined by Bullmore and Sporns (2009), in the process aimed at identifying functional networks, there are many choices and assumptions that can influence the final results. Thus, attention must be paid to the prior knowledge of the system and to the scientific questions we are interested in. If we look at in vitro studies at cellular level, we find that most of these studies are restricted to very small neuronal circuits, typically investigated with patch-clamp techniques (Bi and Poo, 1999), or to large complex networks in which (functional or effective) connectivity cannot be reliably estimated based on statistical measurements (Garofalo et al., 2009) from extracellular recordings. This is particularly true when such electrophysiological recordings imply a large spatial under-sampling, as it is the case with conventional multi-electrode arrays (Maccione et al., 2010). Furthermore, no other experimental investigations or theoretical considerations can be easily made about how networks connect when scaling from local to global interaction and about how information is processed at different levels of complexity (i.e., from the firing of a few interconnected cells up to a coordinated-orchestrated and collective synchronized behavior). In our approach, in an attempt to overcome some of the above issues, we used low-density cultures and we recorded their spontaneous activity by exploiting the high-spatial and temporal resolution allowed by the high-density MEAs. In such experimental configuration we could perform a direct estimation of functionaleffective relationships by using simple CC-based methods. The implemented methods do not allow to directly estimate causal relationships but they consistently preserve those connections that could be sustained by causal interactions. We could identify small directed graphs (Figs. 5 and 6) and we could analyze some of the connectivity features and peculiarities of two and three interconnected neurons (Fig. 7). We extracted such functional-effective interconnections by applying a CC based algorithm and additional filtering procedures, on spontaneous electrophysiological recordings and we presented the estimated maps superimposed to the immuno-fluorescent images of the network structure. Based on previously validated methods, on simulated data (Garofalo et al., 2009), we first showed that such functional maps are stable and that a recording duration of 10 min, even with low firing regimes,

is enough to give a good statistical set to estimate connections. Second, we analyzed the characteristics of these estimated maps in terms of connections strength, delay and length. Third, by taking advantage of the high spatial resolution recordings, we scaled down from the global network connectivity to the identification of simple microcircuits constituted by a few or several interconnected neurons. Those microcircuits, shown superimposed to the fluorescence images, let us estimate more specific and local interactions constituted by two or three neurons. It should be also underlined that in the implementation of the statistical methods we took a restrictive and conservative approach (cf., Section 2.4), by focusing primarily on the main pathways and connections among the statistically plausible ones (i.e., only the strongest connections are considered). Furthermore, we considered only connections compatible with physiological constraints (i.e., propagation velocity) and then we tried to separate efficient and fast mono synaptic links from less reliable and possibly poly-synaptic links (Figs. 5 and 6). In a particular experiment (hn2), a cluster of cell bodies is clearly visible in the bottom part of the culture. This unwanted situation, allowed us to test the capability of the method to individuate tight and short range interactions among cells. In the close-up presented in Fig. 6B, many unresolved bi-directional links are visible, demonstrating highly synchronous coordinated activity. We cannot exclude neither electrically mediated synaptic interactions (i.e., gap junction) nor chemical synapses (many links also with peak delays greater than 1 or 2 ms). Nevertheless, the CC based method well defines such cluster with synchronous activity in terms of functionally hyper-connected graphs. With increased temporal resolution (available with APS technology) it would have been possible try to discern also those propagation directions. This could have been afforded by zooming on a region of interest when recording with the APS-MEA but it would have required having an a priori knowledge of the region of interest. Successively, by scaling down from the entire network to localized microcircuits we considered first neuron’s pair and then propagating pathways constituted by three neurons forming a closed loop. For neuron’s pair analysis, we fit the CC with a Gaussian in order to verify that the CC selected by the STF 0–1 and by the STF greater than 1 ms were characterized by a comparable mean variance. This demonstrates that actual synaptic effective connectivity is found and beside the expected positive correlation between transmission delays and mean variance values, it proves also that fast transmission is efficiently sustained by integrating processes at the basis of synaptic connectivity. The reason we think we can exclude we are measuring intra-neuronal propagation is given by the nondeterministic responses observed in the supposed post-synaptic neurons and by the clear identification of cell bodies in correspondence of the recording site by the immunofluorescence images (e.g., Fig. 7B). For microcircuits made up of three interconnected neurons, we found that, even if we did not apply any partialization on the CC functions (Eichler et al., 2003), the vast majority of the direct connection delays (between neuron 1 and 3) are shorter than the indirect ones (between neuron 1 and 3 via neuron 2) path (cf., Fig. 7D), thus showing that, after the application of the PF filter, actual synaptic connections have been correctly identified in such microcircuits. We think that the proposed approach could support new investigations along the direction of a better understanding of the interplay between structure and function, filling the gap from single neuron level to complex and large network level. The study is still at the beginning and more quantitative evaluations, especially from the morphological–structural point of view are needed. This is certainly beyond the scope of the presented work but the use of patterned networks would definitely allow to design

A. Maccione et al. / Journal of Neuroscience Methods 207 (2012) 161–171

specific microcircuits (with specific topology) and to study, with single cell definition the relationships between structure, topology and function and, as a consequence, to better address the issue of how information is encoded by those microcircuits. References Berdondini L, Imfeld K, Maccione A, Tedesco M, Neukom S, Koudelka-Hep M, et al. Active pixel sensor array for high spatio-temporal resolution electrophsyiological recordings from single cell to large scale neuronal networks. Lab Chip 2009a;9:2644–51. Berdondini L, Massobrio P, Chiappalone M, Tedesco M, Imfeld K, Maccione A, et al. Extracellular recordings from locally dense microelectrode arrays coupled to dissociated cortical cultures. J Neurosci Methods 2009b;177:386–96. Bettencourt LMA, Stephens GJ, Ham MI, Gross GW. Functional structure of cortical neuronal networks grown in vitro. Phys Rev 2007:12. Bi G, Poo M. Synaptic modification by correlated activity: Hebb’s postulate revisited. Annu Rev Neurosci 2001;24:139–66. Bi GQ, Poo MM. Distributed synaptic modification in neural networks induced by patterned stimulation. Nature 1999;401:792–6. Bonifazi P, Ruaro ME, Torre V. Statistical properties of information processing in neuronal networks. Eur J Neurosci 2005;22:2953–64. Brewer GJ, Torricelli JR, Evege EK, Price PJ. Optimized survival of hippocampal neurons in B27-supplemented Neurobasal, a new serum-free medium combination. J Neurosci Res 1993;35:567–76. Bullmore E, Sporns O. Complex brain networks: graph theoretical analysis of structural and functional systems. Nature 2009;10:13. Chiappalone M, Novellino A, Vajda I, Vato A, Martinoia S, van Pelt J. Burst detection algorithms for the analysis of spatio-temporal patterns in cortical networks of neurons. Neurocomputing 2005;65–66:653–62. Duda RO, Hart PE, Stork DG. Pattern classification. 2nd ed. Wiley; 2000. Eichler M, Dahlhaus R, Sandkuhler J. Partial correlation analysis for the identification of synaptic connections. Biol Cybern 2003;89:289–302. Eytan D, Minerbi A, Ziv N, Marom S. Dopamine-induced dispersion of correlations between action potentials in networks of cortical neurons. J Neurophysiol 2004;92:1817–24. Feldt S, Bonifazi P, Cossart R. Dissecting functional connectivity of neuronal microcircuits: experimental and theoretical insight. Trends Neurosci 2011;34:225–36. Garofalo M, Nieus T, Massobrio P, Martinoia S. Evaluation of the performance of information theory-based methods and cross-correlation to estimate the functional connectivity in cortical networks. PLoS One 2009;4:e6482.

171

Grun S, Rotter S. Analysis of parallel spike trains. Springer Series in Computational Neuroscience, vol. 7; 2010. p. 443. Harris KD. Neural signatures of cell assembly organization. Nat Rev Neurosci 2005;6:399–407. Imfeld K, Garenne A, Neukom S, Maccione A, Martinoia S, Koudelka-Hep M, et al. High-resolution MEA platform for in vitro electrogenic cell networks imaging. Conf Proc IEEE Eng Med Biol Soc 2007;2007:6086–9. Imfeld K, Neukom S, Maccione A, Bornat Y, Martinoia S, Farine P-A, et al. High-resolution data acquisition system for extracellular recording of electrophysiological activity. IEEE Trans Biomed Eng 2008:55. Knox CK. Detection of neuronal interactions using correlation analysis. Trends Neurosci 1981;4:222–5. Maccione A, Gandolfo M, Massobrio P, Novellino A, Martinoia S, Chiappalone M. A novel algorithm for precise identification of spikes in extracellularly recorded neuronal signals. J Neurosci Methods 2009;177:241–9. Maccione A, Gandolfo M, Tedesco M, Nieus T, Imfeld K, Martinoia S, et al. Experimental investigation on spontaneously active hippocampal cultures recorded by means of high-density MEAs: analysis of the spatial resolution effects. Front Neuroeng 2010:3. Nicolelis MAL, Lebedev MA. Principles of neural ensemble physiology underlying the operation of brain–machine interfaces. Nat Rev Neurosci 2009;10: 530–40. Pajevic S, Plenz D. Efficient network reconstruction from dynamical cascades identifies small-world topology of neuronal avalanches. PLoS Comput Biol 2009;5:e1000271. Rieke F, Warland D, de Ruyter van Steveninck R, Bialek W. Spikes: exploring the neural code. Cambridge, MA: The MIT Press; 1997. Salinas E, Sejnowski TJ. Correlated neuronal activity and the flow of neural information. Nat Rev Neurosci 2001;2:539–50. Shefi O, Golding I, Segev R, Ben-Jacob E, Ayali A. Morphological characterization of in vitro neuronal networks. Phys Rev E Stat Nonlin Soft Matter Phys 2002;66:021905. Sporns O. Small-world connectivity, motif composition, and complexity of fractal neuronal connections. Biosystems 2006;85:55–64. Sporns O, Tononi G. Classes of network connetivity and dynamics. Complexity 2002;7:28–38. Sporns O, Tononi G, Kötter R. The human connectome: a structural description of the human brain. PLoS Comput Biol 2005;1(4)., http://dx.doi.org/10.1371/journal.pcbi.0010042. Takahashi N, Sasaki T, Matsumoto W, Matsuki N, Ikegaya Y. Circuit topology for synchronizing neurons in spontaneously active networks. Proc Natl Acad Sci USA 2010;107:10244–9.