A systematic assessment of global neural network connectivity through direct electrophysiological assays has remained technically infeasible, even in simpler systems like dissociated neuronal cultures. We introduce an improved algorithmic approach based on Transfer Entropy to reconstruct structural connectivity from network activity monitored through calcium imaging. We focus in this study on the inference of excitatory synaptic links. Based on information theory, our method requires no prior assumptions on the statistics of neuronal firing and neuronal connections. The performance of our algorithm is benchmarked on surrogate time series of calcium fluorescence generated by the simulated dynamics of a network with known ground-truth topology. We find that the functional network topology revealed by Transfer Entropy depends qualitatively on the time-dependent dynamic state of the network (bursting or non-bursting). Thus by conditioning with respect to the global mean activity, we improve the performance of our method. This allows us to focus the analysis to specific dynamical regimes of the network in which the inferred functional connectivity is shaped by monosynaptic excitatory connections, rather than by collective synchrony. Our method can discriminate between actual causal influences between neurons and spurious non-causal correlations due to light scattering artifacts, which inherently affect the quality of fluorescence imaging. Compared to other reconstruction strategies such as cross-correlation or Granger Causality methods, our method based on improved Transfer Entropy is remarkably more accurate. In particular, it provides a good estimation of the excitatory network clustering coefficient, allowing for discrimination between weakly and strongly clustered topologies. Finally, we demonstrate the applicability of our method to analyses of real recordings of in vitro disinhibited cortical cultures where we suggest that excitatory connections are characterized by an elevated level of clustering compared to a random graph (although not extreme) and can be markedly non-local.
The identification of the topological features of neuronal circuits is an essential step towards understanding neuronal computation and function. Despite considerable progress in neuroanatomy, electrophysiology and imaging [1]–[8], the detailed mapping of neuronal circuits is already a difficult task for a small population of neurons, and becomes impractical when accessing large neuronal ensembles. Even in the case of cultures of dissociated neurons, in which neuronal connections develop de novo during the formation and maturation of the network, very few details are known about the statistical features of this connectivity, which might reflect signatures of self-organized critical activity [9]–[11]. Neuronal cultures have emerged in recent years as simple, yet versatile model systems [12], [13] in the quest for uncovering neuronal connectivity [14], [15] and dynamics [16]–[19]. The fact that relatively simple cultures already exhibit a rich repertoire of spontaneous activity [18], [20] make them particularly appealing for studying the interplay between activity and connectivity. The activity of hundreds to thousands of cells in in vitro cultured neuronal networks can be simultaneously monitored using calcium fluorescence imaging techniques [14], [21], [22]. Calcium imaging can be applied both in vitro and in vivo and has the potential to be combined with stimulation techniques like optogenetics [23]. A major drawback of this technique, however, is that the typical frame rate during acquisition is slower than the cell' s firing dynamics by an order of magnitude. Furthermore the poor signal-to-noise ratio makes the detection of elementary firing events difficult. Neuronal cultures are unique platforms to investigate and quantify the accuracy of network reconstruction from activity data, extending analysis tools initially devised for the characterization of macro-scale functional networks [24], [25] to the micro-scale of a developing local circuit. Here we report a new technique based on information theory to reconstruct the connectivity of a neuronal network from calcium imaging data. We use an extension of Transfer Entropy (TE) [26]–[28] to extract a directed functional connectivity network in which the presence of a directed edge between two nodes reflects a direct causal influence by the source to the target node [29]–[31]. Note that “causal influence” is defined operationally as “improved predictability” [32], [33] reflecting the fact that knowledge of the activity of one node (putatively pre-synaptic) is helpful in predicting the future behavior of another node (putatively post-synaptic). TE has previously been used to study gene regulatory networks [34], the flow of information between auditory neurons [35], to infer directed interactions between brain areas based on EEG recordings [36] or between different LFP frequency bands [37], as well as for the reconstruction of the connectivity based on spike times [38], [39]. Importantly, our data-driven TE approach is model-independent. This is in contrast with previous approaches to network reconstruction, which were most often based on the knowledge of precise spike times [40]–[45], or explicitly assumed a specific model of neuronal activity [43], [44]. A problem inherent to the indirect algorithmic inference of network connectivity from real data is that the true target topology of the network is not known and that, therefore, it is difficult to assess the quality of the reconstruction. In order to characterize the behavior of our algorithm and to benchmark its potential performance, we resort therefore to synthetic calcium fluorescence time series generated by a simulated cultured neural network that exhibits realistic dynamics. Since the “ground truth” topology of cultures in silico is known and arbitrarily selectable, the quality of our reconstruction can be evaluated by systematically comparing the inferred with the real network connectivities. We use a simplified network simulation to generate surrogate imaging data, improving their realism with the reproduction of light scattering artifacts [6] which ordinarily affect the quality of the recording. Our surrogate data also reproduce another general feature of the activity of neuronal cultures, namely the occurrence of temporally irregular switching between states of asynchronous activity, with relatively weak average firing rates, and states of highly synchronous activity, commonly denoted as “network bursts” [20], [46], [47]. This switching dynamics poses potentially a major obstacle to reconstruction, since directed functional connectivity can be very different during bursting and inter-burst phases and can bear a resemblance to the underlying structural (i. e. synaptic) connectivity only in selected dynamical regimes in which causal influences reflect dominantly mono-synaptic interactions. To restrict our analysis to such “good” regimes, we resort to conditioning with respect to the averaged fluorescence level, as an indirect but reliable indicator of the network collective dynamics. Appropriate conditioning —combined with a simple correction coping with the poor time-resolution of imaging data— allows the method to achieve a good topology reconstruction performance (assessed from synthetic data), out-performing other standard approaches, without the need to infer exact spike times through sophisticated techniques (as is required, on the contrary, in [43], [45]). Finally, we apply our algorithm —optimized through model-based validation— to the analysis of real calcium imaging recordings. For this purpose, we study spontaneously developing networks of dissociated cortical neurons in vitro and we address, as a first step toward a full topology reconstruction, the simpler problem of extracting only their excitatory connectivity. Early mature cultures display a bursting dynamics very similar to our simulated networks, with which they also share an analogous state-dependency of directed functional connectivity. Our generalized TE approach thus identifies network topologies with characteristic and non-trivial features, like the existence of non-local connections, a broadened and strongly right-skewed distribution of degrees (although not “scale free”) and a moderate but significant level of clustering.
In this study, we consider recordings from in vitro cultures of dissociated cortical neurons (see Materials and Methods). To illustrate the quality of our recordings, in Figure 1A we provide a bright field image of a region of a culture together with its associated calcium fluorescence. As previously anticipated, to simplify the network reconstruction problem, experiments are carried out with blocked inhibitory GABA-ergic transmission, so that the network activity is driven solely by excitatory connections. We record activity of early mature cultures at day in vitro (DIV) 9–12. Such young but sufficiently mature cultures display rich bursting events, combined with sparse irregular firing activity during inter-burst periods (cfr. Discussion). In Figure 1B (left panel) we show actual recordings of the fluorescence traces associated to five different neurons. The corresponding population average for the same time window is shown in Figure 1C (left panel). In these recordings, a stable baseline is broken by intermittent activity peaks that correspond to synchronized network bursts recruiting many neurons. The bursts display a fast rise of fluorescence at their onset followed by a slow decay. In addition, during inter-burst periods, smaller modulations above the baseline are sometimes visible despite the poor time-resolution of a frame every few tens of milliseconds. In order to benchmark and optimize different reconstruction methods, we also generate surrogate calcium fluorescence data (shown in the right panels of Figures 1B and 1C), based on the activity of simulated networks whose ground truth topology is known. We simulate the spontaneous spiking dynamics of networks formed by excitatory integrate-and-fire neurons, along a duration of 60 minutes of real time, matching typical lengths of actual recordings. Calcium fluorescence time series are then produced based on this spiking dynamics, resorting to a model introduced in [43] and described in the Materials and Methods section. Although over 1000 cells are accessible in our experiments, we observed in the simulations that neurons suffice to reproduce the same dynamical behavior observed for larger network sizes, while still allowing for an exhaustive exploration of the entire algorithmic parameter space. Furthermore, despite their reduced density, we maintain in our simulated cultures the same average probability of connection as in actual cultures, where this probability (, see Materials and Methods) is an estimate based on independent studies [15], [48]. The fluorescence signal of a particular simulation run or experiment can be conveniently studied in terms of the distribution of fluorescence amplitudes. As shown in Figure 1D for both simulations and experiments, the amplitude distributions display a characteristic right-skewed shape that emerge from the switching between two distinct dynamical regimes, namely the presence or absence of bursts. The distribution in the low fluorescence region assumes a Gaussian-like shape, corresponding to noise-dominated baseline activity, while the high fluorescence region displays a long tail with a cut-off at the level of calcium fluorescence of the highest network spikes. As we will show later, qualitative similarity between the shapes of the simulated and experimental fluorescence distributions will play an important guiding role for an appropriate network reconstruction. Neurons grown in vitro develop on a two-dimensional substrate and, hence, both connectivity and clustering may be strongly sensitive to the physical distance between neurons. At the same time, due to long axonal projections [13], [49], excitatory synaptic connections might be formed at any distance within the whole culture and both activity and signaling-dependent mechanisms might shape non-trivially long-range connectivity [50], [51]. To test the reconstruction performance of our algorithm, we consider two general families of network topologies that cover a wide range of clustering coefficients. In a first one, clustering occurs between randomly positioned nodes (non-local clustering). In a second one, the connection probability between two nodes decays with their Euclidean distance according to a Gaussian distribution and, therefore, connected nodes are also likely to be spatially close. In particular, in this latter case, the overall level of clustering is determined by how fast the connection probability decays with distance (local clustering). Cortical slice studies revealed the existence of both local [52], [53] and non-local [7], [54] types of clustering. We will later benchmark reconstruction performance for both kinds of topologies and for a wide range of clustering levels, because very similar patterns of neuronal activity can be generated by very different networks, as we now show. Figure 2 illustrates the dynamic behavior of three networks (in this case from the non-local clustering ensemble). The networks are designed to have different clustering coefficients but the same total number of links (see the insets of Figure 2B for an illustration). The synaptic coupling between neurons was adjusted in each network using an automated procedure to obtain bursting activities with comparable bursting rates (see Materials and Methods for details and Table 1 for the actual values of the synaptic weight). As a net effect of this procedure, the synaptic coupling between neurons is slightly reduced for larger clustering coefficients. The simulated spiking dynamics is shown in the raster plots of Figure 2A. These three networks display indeed very similar bursting dynamics, not only in terms of the mean bursting rate, but also in terms of the entire inter-burst interval (IBIs) distribution, shown in Figure 2B. In the same manner, we constructed and simulated local networks —with a small length scale corresponding to high clustering coefficients and vice versa— and obtained the same result, i. e. very similar dynamics for very different decay lengths (not shown). We stress that our procedure for the automatic generation of networks with similar bursting dynamics was not guaranteed to converge for such a wide range of clustering coefficients. Thus, the illustrative simulations of Figure 2 provide genuine evidence that the relation between network dynamics and network structural clustering is not trivially “one-to-one”, despite the fact that more clustered networks have been shown to have different cascading dynamics at the onset of a burst [42]. We focus, first, on the reconstruction of simulated networks, taken from the local and non-local ensembles described above. We compute their directed functional connectivity based on simulated calcium signals. Synthetic fluorescence time series are pre-processed only by simple discrete differentiation, such as to extract baseline modulations associated to potential firing. These differentiated signals are then used as input to any further analyses. Immediately prior to the onset of a burst the network is very excitable. In such a situation it is intuitive to consider that the directed functional connectivity can depart radically from the structural excitatory connectivity, because local events can potentially induce changes at very long ranges due to collective synchronization rather than to direct synaptic coupling. Conversely, in the relatively quiet inter-burst phases, a post-synaptic spike is likely to be influenced solely by the presynaptic firing history. Hence, the directed functional connectivity between neurons is intrinsically state dependent (cfr. also [55]), a property that must be taken into account when reconstructing the connectivity. We illustrate here the state dependency of directed functional connectivity by generating a random network from the local clustering ensemble and by simulating its dynamics, including light scattering artifacts to obtain more realistic fluorescence signals. The resulting distribution of fluorescence amplitudes is divided into seven non-overlapping ranges of equal width, each of them identified with a Roman numeral (Figure 3A). Finally, TE is computed separately for each of these ranges, based on different corresponding subsets of data from the simulated recordings. For simulated data, the inferred connectivity can be directly compared to the ground truth, and a standard Receiver-Operator Characteristic (ROC) analysis can be used to quantify the quality of reconstruction. ROC curves are generated by gradually moving a threshold level from the lowest to the highest TE value, and by plotting at each point the fraction of true positives as a function of the fraction of false positives. The quality of reconstruction is then summarized in a single number by the performance level, which, following an arbitrary convention, is measured as the fraction of true positives at 10% of false positives read out of a complete ROC curve. We plot the performance level as a function of the average fluorescence amplitude in each interval, as shown by the blue line of Figure 3B. The highest accuracy is achieved in the lowest fluorescence range, denoted by I, and reaches a remarkably elevated value of approximately 70% of true positives. The performance in the higher ranges II to IV decreases to a value around 45%, to abruptly drop at range V and above to a final plateau that corresponds to the 10% performance of a random reconstruction (ranges VI and VII). Note that fluorescence values are not distributed homogeneously across ranges I–VII, as evidenced by the overall shape of the fluorescence distribution in Figure 3A. For example, the lowest and highest ranges (I and VII) differ by two orders of magnitude in the number of data points. To discriminate unequal-sampling effects from actual state-dependent phenomena, we studied the performance level using an equal number of data points in all ranges. Effectively, we restrict the number of data points available in each range to be equal to the number of samples in the highest range, VII. The quality of such a reconstruction is shown as the red curve in Figure 3B. The performance level is now generally lower, reflecting the reduced number of time points which are included in the analysis. Interestingly, the “true” peak of reconstruction quality is shifted to range II, corresponding to fluorescence levels just above the Gaussian in the histogram of Figure 3A. This range is therefore the most effective in terms of reconstruction performance for a given data sampling. For the ranges higher than II, the reconstruction quality gradually decreases again to the 10% performance of purely random choices in ranges VI and VII. The effect of adopting a (shorter) equal sample size is particularly striking for range I, which drops from the best performance level almost down to the baseline for random reconstruction. As a matter of fact, range I is the one for which the shrinkage of sample length due to the constraint for uniform data sampling is most extreme (see later section on dependence of performance from sample size). The above analysis leads to a different functional network for each dynamical range studied. For the analysis with an equal number of data point per interval, the seven effective networks are drawn in Figure 3C (for clarity only the top 10% of links are shown). Each functional network is accompanied with the corresponding ROC curve. The lowest range I corresponds to a regime in which spiking-related signals are buried in noise. Correspondingly, the associated functional connectivity is practically random, as indicated by a ROC curve close to the diagonal. Nevertheless, information about structural topology is still conveyed in the activity associated to this regime and can be extracted through extensive sampling. At the other extreme, corresponding to the upper ranges V to VII —associated to fully developed synchronous bursts— the functional connectivity has also a poor overlap with the underlying structural network. As addressed later in the Discussion section, functional connectivity in regimes associated to bursting is characterized by the existence of hub nodes with an elevated degree of connection. The spatio-temporal organization of bursting can be described in terms of these functional connectivity hubs, since nodes within the neighborhood of the same functional hub provide the strongest mutual synchronization experienced by an arbitrary pair of nodes across the network (see Discussion and also Figure S2). The best agreement between functional and excitatory structural connectivity is clearly obtained for the inter-bursts regime associated with the middle range II, and to a lesser degree in ranges III and IV, corresponding to the early building-up of synchronous bursts. Overall, this study of state-dependent functional connectivity provides arguments to define the optimal dynamical regime for network reconstruction: The regime should include all data points whose average fluorescence across the population is below a “conditioning level”, located just on the right side of the Gaussian part of the histogram of the average fluorescence (see Materials and Methods). This selection excludes the regimes of highly synchronized activity (ranges III to VII) and keeps most of the data points for the analysis in order to achieve a good signal-to-noise ratio. Thus, the inclusion of both ranges I and II combines the positive effects of correct state selection and of extensive sampling. The state-dependency of functional connectivity is not limited to synthetic data. Very similar patterns of state-dependency are observed also in real data from neuronal cultures. In particular, in both simulated and real cultures, the functional connectivity associated to the development of bursts displays a stronger clustering level than in the inter-burst periods. An analysis of the topological properties of functional networks obtained from real data in different states (compared with synthetic data) is provided in Figure S3. In this same figure, sections of fluorescence time-series associated to different dynamical states are represented in different colors, for a better visualization of the correspondence between states and fluorescence values (for simplicity, only four fluorescence ranges are distinguished). Our generalized TE, conditioned to the proper dynamic range, enables the reconstruction of network topologies even in the presence of light scattering artifacts. For non-locally clustered topologies we obtain a remarkably high accuracy of up to 75% of true positives at a cost of 10% of false positives. An example of the reconstruction for such a network, with, is shown in Figure 4A. For locally-clustered topologies, accuracy typically reaches 60% of true positives at a cost of 10% of false positives, and an example for is shown in Figure 5A. In both topologies, we observe that for a low fraction of false positives detection (i. e. at high thresholds) the ROC curve displays a sharp rise, indicating a very reliable detection of the causally most efficient excitatory connections. A decrease in the slope, and therefore a rise in the detection of false positives and a larger confidence interval, is observed only at higher fractions of false positives. The confidence intervals are broader in the case of locally-clustered topologies because of the additional network-to-network variability that results from the placement of neurons (which is irrelevant for the generation of the non-locally clustered ensembles, see Materials and Methods). The performance level (fraction of true positives for 10% of false positives, denoted by) provides a measure of the quality of the reconstruction, and allows the comparison of different methods for different network topologies, conditioning levels, and external artifacts (i. e. presence or absence of simulated light scattering). We test linear methods, XC and GC (of order 2; the performance of GC of order 1 is very similar and not shown), and non-linear methods, namely MI and TE (of Markov orders 1 and 2; see Materials and Methods for details). XC and MI are correlation measures, while GC and TE are causality measures. Note that, for each of these methods, we account for state dependency of functional connectivity, performing state separation as described in the Materials and Methods section. In the case of the non-local clustering ensemble and without light scattering (Figure 6, top row), even a linear method such as XC achieves a good reconstruction. This success indicates an overlap between communities of higher synchrony in the calcium fluorescence, associated to stronger activity correlations, and the underlying structural connectivity, especially for higher full clustering indices. GC-based reconstructions have an overall worse quality, due to the inadequacy of a linear model for the prediction of our highly nonlinear network dynamics, but they show similarly improved performance for higher. In a band centered around a shared optimal conditioning level, both MI and generalized TE show a robust performance across all clustering indices. This value is similar to the upper bound of the range II depicted in Figure 3A, i. e. it lies at the interface between the bursting and silent dynamical regimes. In particular for TE and in the case of low clustering indices (which leads to networks closer to random graphs), conditioning greatly improves reconstruction quality. At higher clustering indices the decay in performance is only moderate for conditioning levels above the optimal value, indicating an overlap between the functional connectivities in the bursting and silent regimes. Note, on the contrary, that the performance of MI rapidly decreases if a non-optimal conditioning level is assumed. The introduction of light scattering causes a dramatic drop in performance of the two linear methods (XC and GC), and even of MI and TE with Markov order. The performance of TE at Markov order 2 also deteriorates, but is still significantly above the random reconstruction baseline in a broad region of parameters. Interestingly, for the optimal conditioning level the performance of the TE for does not fall below for any clustering level or value. It is precisely in this optimal conditioning range that we obtain the linear relations between reconstructed and structural clustering coefficients, for both the non-local and the local clustering ensembles. A similar trend is obtained when varying the length scale in the local ensembles (see Supplementary Figure S5). For very local clustering and without light scattering, both XC and TE achieve performance levels up to 80%. The introduction of light scattering, however, reduces the performance of all measures except for MI (but only in the narrow optimal conditioning range) and for TE of higher Markov orders (robust against non-optimal selection of conditioning level). Overall, the performance of the reconstruction for the local clustering ensembles is lower than for the non-locally clustered ensembles. This is also true, incidentally, in absence of light scattering since networks sampled from this ensemble tend to be very similar to purely random topologies (of the Erdös-Rényi type, see e. g. [62]) as soon as the length scale is sufficiently long, and for which performance is generally poorer (cfr. top row of Figure 6, for weak clustering levels). Our new TE method significantly improves the reconstruction performance compared to the original TE formulation [Eq. (9) ]. As shown in Figure 7A for both the local and the non-local clustered networks, reconstruction with the original TE formulation (Figure 7A, blue line) yields worse results than a random reconstruction, as indicated by the corresponding ROC curves falling below the diagonal. Such a poor performance is due in large part to “misinterpreted” delayed interactions. Indeed, by taking into account same bin interactions, a boost in performance is observed (red line). Figure 7A also shows that an additional leap in performance is obtained when the analysis is conditioned (i. e. restricted) to a particular dynamical state of the network, increasing reconstruction quality by 20% (yellow line in Figure 7A). The determination of the optimal conditioning level is discussed later and takes into account the considerations introduced above (cfr. Figure 3). Note that the introduction of same bin interactions alone (red color curves) or conditioning on the dynamical state of network alone (yellow color curves) already brings the performance to a level well superior to random performance. However, at least for our simulated calcium-fluorescence time series, a remarkable boost in performance is obtained only when the inclusion of same-bin interactions and optimal conditioning are combined together (green color curves). Although, in principle, conditioning is enough to indirectly select a proper dynamical regime, the poor time-resolution of the analyzed signals (constrained not only by the frame-rate of acquisition but also intrinsically by the kinetics of the dissociation reaction of the calcium-sensitive dye [63]) also requires the potential consideration of causally-linked events occurring in the same time-bin. A different way to represent reconstruction performance are “Positive Precision Curves” (as introduced in [56] and described in the Materials and Methods section), obtained by plotting, at a given number of reconstructed links, the “true-false ratio” (TPR), which emphasizes the probability that a reconstructed link is present in the ground truth topology (true positive). For the same networks and reconstruction as above, we plot the PPCs in Figure S6 (for the ROC curves see Figure 7A). Over a wide range of the number of reconstructed links (TFS), the PPC displays positive values of the TFR, indicative of a majority of true positives over false positives. For both the locally and non-locally clustered networks, the PPC reaches a maximum value of the TPR about 0. 5 and remains positive up to about 18% of included links for the clustered topology, or up to 12% in the case of the local topology. In Figure 7B, we analyze the performance of our algorithm against changes of the sample size. Starting from simulated recordings lasting 1 h of real time (corresponding to about 360 bursting events) and with a full sample number of, we trimmed these recordings producing shorter fluorescence time series with samples, with being a divisor of the sample size. For both network topology ensembles, we found that a reduction in the number of samples by a factor of two (corresponding to 30 minutes or about 180 bursts) still yields a performance level of. By further reducing the sample size, we reach a plateau with a quality of for about 40 bursts (corresponding to 6 minutes). All the experiments analyzed in this work are carried out with a duration between 30 and 60 minutes. Since conditioning, needed to achieve high performance, requires one to ignore a conspicuous fraction of the recorded data, we expect long recordings to be necessary for a good reconstruction, albeit the fact that it is possible to increase the signal-to-noise ratio by increasing the intensity of the fluorescent light. However, the latter manipulation has negative implications for the health of neurons due to photo-damage, limiting our experimental recordings to a maximum of 2 hours. We apply our analysis to actual recordings from in vitro networks derived from cortical neurons (see Materials and Methods). To simplify the network reconstruction problem, experiments are carried out with blocked inhibitory GABA-ergic transmission, so that the network activity is driven solely by excitatory connections. This is consistent with previously discussed simulations, in which only excitatory neurons were included. We consider in Figure 8 a network reconstruction based on a 60 minutes recording of the activity of a mature culture, at day in vitro (DIV) 12, in which active neurons were simultaneously imaged. A fully analogous network reconstruction for a second, younger dataset at DIV 9 is presented in Supplementary Figure S7. In general, fluorescence data neither affected by photo-bleaching nor by photo-damage during this time, as proved by the stability of the average fluorescence signal shown in the Supplementary Figure S8A. The probability distribution of the average fluorescence signal is computed in the same way as for the simulated data. Neuronal dynamics and the calcium fluorescence display the same bursting dynamics that are well captured by the simulations, leading to a similar fluorescence distribution (Figure 1D). Thanks to this similarity we can make use of the intuition developed for synthetic data to estimate an adequate conditioning level. We select therefore a conditioning level such as to exclude the right-tail of high fluorescence associated to to fully-developed bursting transient regimes. We have verified, however, that the main qualitative topological features of the reconstructed network are left unchanged when varying the conditioning level in a range centered on our “optimal” selection. More details on conditioning level selection are given in the Materials and Methods section. Reconstruction analysis is carried out for the entire population of imaged neurons. We analyze a network defined by the top 5% of TE-ranked links, as discussed in a previous section. Such choice leads to an average in–degree of about 100, compatible with average degrees reported previously for neuronal cultures of corresponding age (DIV) and density [14], [64].
We have introduced a novel extension of Transfer Entropy, an information theoretical measure, and applied it to infer excitatory connectivity of neuronal cultures in vitro. Other studies have previously applied TE (or a generalization of TE) to the reconstruction of the topology of cultured networks [39], [56]. However, our study introduces and discusses important novel aspects, relevant for applications. As shown in Figure 3 for simulated data and in Supplementary Figure S3 for actual culture data, epochs of synchronous bursting are associated to functional connectivity with a stronger degree of clustering and a weaker overlap with the underlying structural topology. This feature of functional connectivity is tightly related to the spatio-temporal organization of the synchronous bursts. In Figure S2A, we highlight the position of selected nodes (of the simulated network considered in Figure 3), characterized by an above-average in-degree of functional connectivity (see Materials and Methods for a detailed definition in terms of TE scores). We denote these nodes —different in general for each of the dynamic regimes numbered from I to VII— as (state-dependent) functional connectivity hubs. Given a specific hub, we can then define the community of its first neighbors in the corresponding functional network. Consistently across all dynamic ranges (but the noise-dominated range I) we find that synchronization within each of these functional communities is significantly stronger than between the communities centered on different hubs ( (, Mann-Whitney test, see Materials and Methods). The results of this comparison are reported in Figure S2B, showing particularly marked excess synchronization during burst build-up (ranges II, III and IV) or just prior to burst waning in the largest-size bursts (VII). Therefore, functional connectivity hubs reflect foci of enhanced local bursting synchrony. Other studies (in brain slices) reported evidence of functional connectivity hubs, whose direct stimulation elicited a strong synchronous activation [69], [70]. In [69], the functional hubs were also structural hubs. In the case of our networks, however only functional hubs associated to ranges II and III have an in-degree (and out-degree) larger than average as in [69]. In the other dynamic ranges, this tight correspondence between structural and functional hubs does not hold anymore. Nevertheless, in all dynamic ranges (but range I), we find that pairs of functional hubs have an approximately three-times larger chance of being structurally connected than pairs of arbitrarily selected nodes (not shown). The timing of firing of these strong-synchrony communities is analyzed in Figure S2C. There we show that the average time of bursting of functional communities for different dynamic ranges is shifted relatively to the average bursting time over the entire network (details of the estimation are provided in Materials and Methods). This temporal shift is negative for the ranges II and III (indicating that functional hubs and related communities fire on average earlier than the rest of the culture) and positive for the ranges V to VII (indicating that firing of these communities occurs on average later than the rest of the culture). The highest negative time delay is detected in range III, such that the communities organized around its associated functional hubs can be described as local burst initiation cores [71], [72]. In this study we did not consider inhibitory interactions, neither in simulations nor in experiments (GABAergic transmission was blocked), but we attempted uniquely the reconstruction of excitatory connectivity. We would like to point out that this is not a general limitation of TE, since the applicability of TE does not rely on assumptions as to the specific nature of a given causal relationship – for instance about whether a synapse is excitatory or inhibitory. In this sense, TE can be seen as a measure for the absolute strength of a causal interaction, and is able in principle to capture the effects on dynamics of both inhibitory and excitatory connections. Note indeed that previous studies [39], [56] have used TE to infer as well the presence of inhibitory connections. However, TE alone could not discriminate the sign of the interaction and additional post-hoc considerations had to be made in order to separate the retrieved connections into separate excitatory and inhibitory subgroups (e. g. , in [39], based on the supposedly known existence of a difference in relative strength between excitatory and inhibitory conductances). In summary, we have developed a new generalization of Transfer Entropy for inferring connectivity in neuronal networks based on fluorescence calcium imaging data. Our new formalism goes beyond previous approaches by introducing two key ingredients, namely the inclusion of same bin interactions and the separation of dynamical states through conditioning of the fluorescence signal. We have thoroughly tested our formalism in a number of simulated neuronal architectures, and later applied it to extract topological features of real, cultured cortical neurons. We expect that, in the future, algorithmic approaches to network reconstruction, and in particular our own method, will play a pivotal role in unravelling not only topological features of neuronal circuits, but also in providing a better understanding of the circuitry underlying neuronal function. These theoretical and numerical tools may well work side by side with new state-of-the-art techniques (such as optogenetics or high-speed two-photon imaging [96]–[100]) that will enable direct large-scale reconstructions of living neuronal networks. Our Transfer Entropy formalism is highly versatile and could be applied to the analysis of in vivo voltage-sensitive dye recordings with virtually no modifications. On a shorter time-scale, it would be important to extend our analysis to the reconstruction of both excitatory and inhibitory connectivity in in vitro cultures, which is technically feasible, and to compare diverse network characteristics, such as neuronal density or aggregation. Our algorithm could be used to systematically reconstruct the connectivity of cultures at different development stages in the quest for understanding the switch from local to global neuronal dynamics. Another crucial open issue is to design suitable experimental protocols allowing to confirm the existence of at least some of the inferred synaptic links, in order to validate statistically the reconstructed connectivity. For instance, the actual presence of directed links to which our algorithm assigns the largest TE scores might be systematically probed through targeted paired electrophysiological stimulation and recording. Furthermore, GFP transfection or inmunostaining might be used to obtain actual, precise anatomical data on network architecture to be compared with the reconstructed one. Finally, it might be interesting to reconstruct connectivity of cultured networks before and after physical disconnection of different areas of the culture (e. g. by mechanical etching of the substrate or by chemical silencing). These manipulations would provide a scenario to verify whether TE-based reconstructions correctly capture the absence of direct connections between areas of the neuronal network which are known to be artificially segregated.
We generated synthetic networks with neurons, distributed randomly over a squared area of 0. 5 mm lateral size. We chose as the connection probability between neurons [48], leading to sparse connectivities similar to those observed in local cortical circuits [53]. We used non-periodic boundary conditions to reproduce eventual “edge” effects that arise from the anisotropic cell density at the boundaries of the culture. We considered two general types of networks: (i) a locally-clustered ensemble, where the probability of connection depended on the spatial distance between two neurons; and (ii) a non-locally clustered ensemble, with the connections engineered to display a certain degree of clustering. For the case of a non-local clustering ensemble, we first created a sparse connectivity matrix, randomly generating links with a homogeneous probability of connection across pairs of neurons. We next selected a random pairs of links and “crossed” them (links and became and). We accepted only those changes that updated the clustering index in the direction of a desired target value, thereby maintaining the number of incoming as well as outgoing connections of each neuron. The crossing process was iterated until a clustering index higher or equal to the target value was reached. The overall procedure led to a full clustering index of the reference random network of (mean and standard deviation, respectively, across 6 networks). After the rewiring iterations, we then achieved standard deviations from the desired target clustering value smaller than 0. 1% for all higher clustering indices. We measured the full clustering index of our directed networks according to a common definition introduced by [65]: (1) The binary adjacency matrix is denoted by, with for a link, and zero otherwise. The adjacent matrix provides a complete description of the network topological properties. For instance, the in-degree of a node can be computed as, and the out-degree as. The total number of links of a node is given by the sum of its in-degree and its out-degree (). The number of bidirectional links of a given node (i. e. links between and so that and are reciprocally connected by directed connections) is given by. The adjacency matrix did not contain diagonal entries. Such entries would correspond “toautaptic” links that connect a neuron with itself. Note that our directed functional connectivity analysis is based on bivariate time series, and therefore it would be structurally unfit to detect this type of links. For the case of the local clustering ensemble, two neurons separated a Euclidean distance were randomly connected with a distance dependent probability described by a Gaussian distribution, of the form, with a characteristic length scale. To guarantee that a constant average number of links was present in the network, this Gaussian distribution was rescaled by a constant pre-factor, obtained as follows. We first generated a network based on the unscaled kernel and computed the resulting number of links. With this value we then generated a final network based on the rescaled kernel. The dynamics of the generated neuronal networks was studied using the NEST simulator [101], [102]. We modeled the neurons as leaky integrate-and-fire neurons, with the membrane potential of a neuron described by [103], [104]: (2) where is the leak conductance and is the membrane time-constant. The term accounts for a time-dependent input current that arises from recurrent synaptic connections. In the absence of synaptic inputs, the membrane potential relaxes exponentially to a resting level set arbitrarily to zero. Stimulation in the form of inputs from other neurons increase the membrane potential, and above the threshold an action potential is elicited (neuronal firing). The membrane voltage is then reset to zero for a refractory period of. The generated action potential excites post-synaptic target neurons. The total synaptic currents are then described by (3) where is the adjacency matrix, and is a synaptic time constant. The resulting excitatory post-synaptic potentials (EPSPs) have a standard difference-of-exponentials time-course [105]. Neurons in culture show a rich spontaneous activity that originates from both fluctuations in the membrane potential and small currents in the pre-synaptic terminals (minis). The latter is the most important source of noise and plays a pivotal role in the generation and maintenance of spontaneous activity [106]. To introduce the spontaneous firing of neurons in Eq. (3), each neuron was driven, through a static coupling conductance with strength, by independent Poisson spike trains (with a stationary firing rate of, spikes fired at stochastic times). Neurons were connected via synapses with short-term depression, due to the finite amount of synaptic resources [104]. We considered only purely excitatory networks to mimic the experimental conditions in which inhibitory transmission is fully blocked. Concerning the recurrent input to neuron, the set represents times of spikes emitted by a presynaptic neuron, is a conduction delay of, while sets a homogeneous scale for the synaptic weights of recurrent connections, whose time-dependent strength depends on network firing history through the equations (4) (5) In these equations, represents the fraction of neurotransmitters in the “effective state”, in the “recovered state” and in the “inactive state” [103], [104]. Once a pre-synaptic action potential is elicited, a fraction of the neurotransmitters in the recovered state enters the effective state, which is proportional to the synaptic current. This fraction decays exponentially towards the inactive state with a time scale, from which it recovers with a time scale. Hence, repeated firing of the presynaptic cell in an interval shorter than gradually reduces the amplitude of the evoked EPSPs as the synapse is experiencing fatigue effects (depression). Random networks of integrate-and-fire neurons coupled by depressing synapses are well-known to naturally generate synchronous events [104], comparable to the all-or-none behavior that is observed in cultured neurons [46], [47]. To obtain in our model a realistic bursting rate [47], the synaptic weight of internal connections was set to result into a network bursting of for all the network realizations that we studied, and in particular for any considered (local or non-local) clustering level. Therefore, after having generated each network topology, we assigned the arbitrary initial value of to internal synaptic weights and simulated 200 seconds of network dynamics, evaluating the resulting average bursting rate. If it was larger (smaller) than the target bursting rate, then the synaptic weight was reduced (increased) by 10%. We then iteratively adjusted by (linearly) extrapolating the last two simulation results towards the target bursting rate, until the result was closer than 0. 01 Hz to the target value. The resulting used values of are provided in Table 1. Note that we defined a network burst to occur when more than 40% of the neurons in the network were active within a time window of 50 ms. Such a criterion does not play any role in the reconstruction algorithm itself, where state selection is achieved through conditioning, but is only used for the automated generation of random networks with a prescribed bursting rate. Typically, for a fully developed burst, more than 90% of the neurons fire within a 50 ms bin, while, during inter-burst intervals, less than 10% do. Due to the clear separation between these two regimes, our burst detection procedure does not depend significantly on the precise choice of threshold within a broad interval. To reproduce the fluorescence signal measured experimentally, we treated the simulated spiking dynamics to generate surrogate calcium fluorescence signals. We used a common model introduced in [43] that gives rise to an initial fast increase of fluorescence after activation, followed by a slow decay (). Such a model describes the intra-cellular concentration of calcium that is bound to the fluorescent probe. The concentration changes rapidly by a step amount of for each action potential that the cell is eliciting in a time step, of the form (6) where is the total number of action potentials. The net fluorescence level associated to the activity of a neuron is finally obtained by further feeding the Calcium concentration into a saturating static non-linearity, and by adding a Gaussian distributed noise with zero mean: (7) For the simulations, we used a saturation concentration of and noise with a standard deviation of 0. 03. We considered the light scattered in a simulated region of interest (ROI) from surrounding cells. Denoting as the distance between two neurons and and by the scattering length scale (determined by the typical light deflection in the medium and the optical apparatus), the resulting fluorescence amplitude of a given neuron is given by (8) A sketch illustrating the radius of influence of the light scattering phenomenon is given in Figure S9. The scaling factor sets the overall strength of the simulated scattering artifact. Note that light scattered, according to the equation shown above, could be completely corrected using a standard deconvolution algorithm, at least for very large fields of view and a scattering length known with sufficient accuracy. In a real setup however, the relatively small fields of view (on the order of), the inaccuracies in inferring the scattering radius, as well as the inhomogeneities in the medium and on the optical system, make perfect deconvolution not possible. Therefore, artifacts due to light scattering cannot be completely eliminated [107], [108]. The scaling factor, that we arbitrarily assumed to be small and with value, can be seen as a measure of this residual artifact component. In its original formulation [26], for two discrete Markov processes and (here shown for equal Markov order), the Transfer Entropy (TE) from to was defined as: (9) where is a discrete time index and is a vector of length whose entries are the samples of at the time steps, , …, . The sum goes over all possible values of, and. TE can be seen as the distance in probability space (known as the Kullback-Leibler divergence [109]) between the “single node” transition matrix and the “two nodes” transition matrix. As expected from a distance measure, TE is zero if and only if the two transition matrices are identical, i. e. if transitions of do not depend statistically on past values of, and is greater than zero otherwise, signaling dependence of the transition dynamics of on. We use TE to evaluate the directed functional connectivity between different network nodes. In a pre-processing step, we apply a basic discrete differentiation operator to calcium fluorescence time series, as a rather crude way to isolate potential spike events. Thus, given a network node, we define. This pre-processing step also improves the signal-to-noise ratio, thus allowing for a better sampling of probability distributions with a limited number of data points. To adapt TE to our particular problem we need to take into account the general characteristics of the system. We therefore modified TE in two crucial aspects: (10) We then include all data points at time instants in which this average fluorescence is below a predefined threshold parameter, i. e. we consider only the time points that fulfill. We only make an exception that corresponds to the simulations of Figure 3 and Figure S2, where we considered time points that fall within an interval bounded by a higher and a lower thresholds, i. e. . Using these two novel aspects, we have extended the original description of Transfer Entropy [Eq. (9) ] to the following form (11) Probability distributions have to be evaluated as discrete histograms. Hence, the continuous range of fluorescence values (see e. g. the bottom panels of Figure 1) is quantized into a finite number of discrete levels. We typically used a small, a value that we justify based on the observation that the resulting bin width is close to twice the standard deviation of the signal. The presence of large fluctuations, most likely associated to spiking events, is then still captured by such a coarse, almost non-parametric description of fluorescence levels. Generalized TE values are obtained for every possible directed pair of network nodes, and using a fixed threshold level. The set of TE scores are then ranked in ascending order and scaled to fall in the unit range. A threshold is then applied to the rescaled data, so that only those links with scores above are retained in the reconstructed network. A standard Receiver-Operator Characteristic (ROC) analysis is used to assess the quality of the reconstruction by evaluating the number of true positives (reconstructed links that are present in the actual network) or false positives (not present), and for different threshold values [110]. The highest threshold value leads to zero reconstructed links and therefore zero true positives and false positives. At the other extreme, the lowest threshold provides both 100% of true positives and false positives. Intermediate thresholds give rise to a smooth curve of true/false positives as a function of the threshold. The performance of the reconstruction is then measured as the degree of deviation of this curve from the diagonal, and that corresponds to a random choice of connections between neurons. To provide a simple method to compare different reconstructions, we arbitrarily use the quantity, defined as the fraction of true positives for a 10% of false positives, as indicator for the quality of the reconstruction. An alternative to the ROC is the Positive Prediction Curve [56], plotting the “true-false ratio” (TFR) against the number of reconstructed links, called “true-false sum” (TFS). The TFR represents the fraction of true positives relative to the false positives. Denoting by #TP the absolute number of true positives and by #FP the number of false positives, TFR is therefore defined in [56] in the following way: (12) The case corresponds to the case that, for any given reconstructed link, it is on average equally likely that it is in fact a true positive rather than a false positive. To gain further insight into the quality of our reconstruction method, we compare reconstructions based on TE with three other reconstruction strategies, namely cross-correlation, mutual information, and Granger causality. Cross-correlation (XC) reconstructions are based on standard Pearson cross-correlation. The score assigned to each potential link is given by the largest cross-correlogram peak for lags between and, of the form (13) In a similar way, the scores for Mutual Information (MI) reconstructions are evaluated as (14) Analogously to TE, the sum goes over all entries of the joint probability matrix. For the reconstruction based on Granger causality (GC) [33] we first model the signal by least-squares fitting of a univariate autoregressive model, obtaining the coefficients and the residual, (15) In a second step, we fit a second bivariate autoregressive model that includes the potential source signal, and determine the residual, (16) Note that in the latter bivariate regression scheme we take into account “same bin” interactions as for Transfer Entropy (index of the second sum starts at). Given, the covariance matrix of the univariate fit in Eq. (15), and, the covariance matrix of the bivariate fit in Eq. (16), GC is then given by the logarithm of the ratio between their traces: (17) GC analyses were performed at an order. Analyses at yielded however fully analogous performance (not shown). We note that the same pre-processing used for TE is also adopted for all the other analyses. The same holds for conditioning on the value of the average fluorescence, which can be applied simply by only including the subset of samples in which. Connectivity in reconstructed networks is often inhomogeneous, and groups of nodes with tighter internal connectivity are sometimes visually apparent (see e. g. reconstructed topologies in Figure 3C). We do not attempt a systematic reconstruction of network communities [111], but we limit ourselves to the detection of “causal sink” nodes [112], which have a larger than average in-degree. We define this property in terms of the sum of TE from all other nodes to one particular node (), choosing the top 20 nodes for each particular network as selected “hub nodes”. We then analyze the dynamics of these selected hub nodes and of their neighbors. Specifically we define as the subgraph spanned by a given hub node and by its first neighbors. We analyze then the cross-correlogram of the average fluorescence of a given group with the average fluorescence of the whole culture: (18) The -notation indicates that we correlate discretely differentiated average fluorescence time series, rather than the average time series themselves. Indeed, cross-correlograms for these differentiated time series are well modeled by a Gaussian functional form, due to the slow change of the averaged fluorescence compared to the sampling rate. Therefore, we fit a Gaussian to the cross-correlogram: (19) determining thus a cross-correlation amplitude, a cross-correlation peak lag and the standard deviation. The cross-correlation peak lag indicates therefore whether nodes in a given local hub neighborhood fire on average earlier or later than other neurons in the network. Relative strength of synchrony within a local hub neighborhood can be analogously evaluated by computing XCs, as defined in Eq. (13), for all the links within and comparing it with peak XCs over the entire network. Primary cultures of cortical neurons were prepared following standard procedures [14], [113]. Cortices were dissected from Sprague-Dawley embryonic rat brains at 19 days of development, and neurons dissociated by mechanical trituration. Neurons were plated onto 13 mm glass cover slips (Marienfeld, Germany) previously coated overnight with 0. 01% Poly-l-lysine (Sigma) to facilitate cell adhesion. Neuronal cultures were incubated at, 95% humidity, and for 5 days in plating medium, consisting of 90% Eagle' s MEM —supplemented with 0. 6% glucose, 1% 100× glutamax (Gibco), and gentamicin (Sigma) — with 5% heat-inactivated horse serum (Invitrogen), 5% heat-inactivated fetal calf serum (Invitrogen), and B27 (Invitrogen). The medium was next switched to changing medium, consisting of of 90% supplemented MEM, 9. 5% heat-inactivated horse serum, and 0. 5% FUDR (5-fluoro-deoxy-uridine) for 3 days to limit glia growth, and thereafter to final medium, consisting of 90% supplemented MEM and 10% heat-inactivated horse serum. The final medium was refreshed every 3 days by replacing the entire culture well volume. Typical neuronal densities (measured at the end of the experiments) ranged between 500 and. Cultures prepared in these conditions develop connections within 24 hours and show spontaneous activity by day in vitro (DIV) 3–4 [14], [18], [20]. GABA switch, the change of GABAergic response from excitatory to inhibitory, occurs at DIV 6–7 [14], [20]. Neuronal activity was studied at day in vitro (DIV) 9–12. Prior to imaging, cultures were incubated for 60 min in pH-stable recording medium in the presence of 0. 4% of the cell-permeant calcium sensitive dye Fluo-4-AM (Invitrogen). Recording solution includes (in mM): , , , , 10 glucose, and 10 HEPES; pH titrated to 7. 4 and osmolarity to with sucrose. The culture was washed off Fluo-4 after incubation and finally placed in a chamber filled with fresh recording medium. The chamber was mounted on a Zeiss inverted microscope equipped with a 5× objective and a 0. 4× optical zoom. Neuronal activity was monitored through high-speed fluorescence imaging using a Hamamatsu Orca Flash 2. 8 CMOS camera attached to the microscope. Images were acquired at a speed of 100 frames/s (i. e. 10 ms between two consecutive frames), which were later converted to a 20 ms resolution using a sliding window average to match the typical temporal resolution of such recordings [13], [43], [45], [90], [114]. The recorded images had a size of pixels with 256 grey-scale levels, and a final spatial resolution of. This settings provided a final field of view of that contained on the order of 2000 neurons. Activity was finally recorded as a long image sequence of 60 minutes in duration. We verified that the fluorescence signal remained stable during the recording, as shown in Supplementary Figure S8A, and we did not observe neither photo-bleaching of the calcium probe nor photo-damage of the neurons. Before the beginning of the experiment, inhibitory synapses were fully blocked with bicuculline, a antagonist, so that activity was solely driven by excitatory neurons. Since cultures were studied after GABA switch, the blockade of inhibition resulted in an increase of the fluorescence amplitude, which facilitated the detection of neuronal firing, as illustrated in Figure S8B. The image sequence was analyzed at the end of the experiment to identify all active neurons, which were marked as regions of interest (ROIs) on the images. The average grey-level on each ROI along the complete sequence finally provided, for each neuron, the fluorescence intensity as a function of time. Each sequence typically contained on the order of a hundred bursts. Examples of recorded fluorescence signal for individual neurons are shown in Figure 1B. The fluorescence data obtained from recordings of neuronal cultures was analyzed following exactly the same procedures used for simulated data (e. g. processed in a pipeline including discrete differentiation, TE or other metrics evaluation, ranking, and final thresholding to maintain the top 10% of connections). Due to the lack of knowledge of the ground-truth topology, optimal conditioning level cannot be known. However, based on the similarity between experimental and simulated distributions of calcium fluorescence, we select a conditioning level such to exclude the high fluorescence transients associated to fully-developed bursting transients while keeping as many data points as possible. Concretely, this is achieved by taking a conditioning level equal to approximately two standard deviations above the mean of a Gaussian fit to the left peak of the fluorescence histogram. Such a level coincides with the point where, when gradually increasing the conditioning level, the reconstructed clustering index reaches a plateau, i. e. matches indicatively the upper limit of range II in Figure 3. To check for robustness of our reconstruction, we generated alternative reconstructions based on different conditioning levels. For the selected conditioning value, and for both the experimental datasets analyzed (Figures 8 and S7), we verified that the inferred topological features, including notably the average clustering coefficient and connection distance, were stable in a range centered on the selected conditioning value and wide as much as approximately two standard deviations of the fluorescence distribution. To identify statistically significant non-random features of the real cultured networks in exam, we compared the reconstructed topology to two randomizations. A first one consisted in a complete randomization that preserved only the total number of connections in the network, but scrambled completely the source and target nodes. The resulting random ensemble of graphs was an Erdös-Rényi ensemble (see, e. g. [62]) in which each possible link exists with a uniform probability of connection, where is the total number of connections in the reference reconstructed network. A second partial randomization preserved the in-degree distributions only, and was implemented by shuffling the entries of each row of the reconstructed adjacency matrix, internally row-by-row. In this way, the out-degrees of each node were preserved. In both randomization processes, we disallowed diagonal entries. For both randomizations we calculated the in-degree, the distance of connections and the full clustering index for each node, leading to distributions of network topology features that could be compared between the reconstructed network and the randomized ensembles, to identify significant deviations from random expectancy.
Unraveling the general organizing principles of connectivity in neural circuits is a crucial step towards understanding brain function. However, even the simpler task of assessing the global excitatory connectivity of a culture in vitro, where neurons form self-organized networks in absence of external stimuli, remains challenging. Neuronal cultures undergo spontaneous switching between episodes of synchronous bursting and quieter inter-burst periods. We introduce here a novel algorithm which aims at inferring the connectivity of neuronal cultures from calcium fluorescence recordings of their network dynamics. To achieve this goal, we develop a suitable generalization of Transfer Entropy, an information-theoretic measure of causal influences between time series. Unlike previous algorithmic approaches to reconstruction, Transfer Entropy is data-driven and does not rely on specific assumptions about neuronal firing statistics or network topology. We generate simulated calcium signals from networks with controlled ground-truth topology and purely excitatory interactions and show that, by restricting the analysis to inter-bursts periods, Transfer Entropy robustly achieves a good reconstruction performance for disparate network connectivities. Finally, we apply our method to real data and find evidence of non-random features in cultured networks, such as the existence of highly connected hub excitatory neurons and of an elevated (but not extreme) level of clustering.
circuit models
connectomics
calcium imaging
developmental neuroscience
computer science
neuroanatomy
neural networks
computational neuroscience
biology
computerized simulations
neuroscience
neural circuit formation
neuroimaging
Model-Free Reconstruction of Excitatory Neuronal Connectivity from Calcium Imaging Signals