Super-selective reconstruction of causal and direct connectivity with application to in-vitro iPSC neuronal networks

TL;DR: A novel mathematically rigorous, model-free method to map effective - direct and causal - connectivity of neuronal networks from multi-electrode array data and is shown to offer important insights into the functional development of in-vitro iPSC-derived neuronal cultures.

Abstract: Despite advancements in the development of cell-based in-vitro neuronal network models, the lack of appropriate computational tools limits their analyses. Methods aimed at deciphering the effective connections between neurons from extracellular spike recordings would increase utility of in-vitro lo- cal neural circuits, especially for studies of human neural development and disease based on induced pluripotent stem cells (hiPSC). Current techniques allow statistical inference of functional couplings in the network but are fundamentally unable to correctly identify indirect and apparent connections between neurons, generating redundant maps with limited ability to model the causal dynamics of the network. In this paper, we describe a novel mathematically rigorous, model-free method to map effective - direct and causal - connectivity of neuronal networks from multi-electrode array data. The inference algorithm uses a combination of statistical and deterministic indicators which, first, enables identification of all existing functional links in the network and then, reconstructs the directed and causal connection diagram via a super-selective rule enabling highly accurate classification of direct, indirect and apparent links. Our method can be generally applied to the functional characterization of any in-vitro neuronal networks. Here, we show that, given its accuracy, it can offer important insights into the functional development of in-vitro iPSC-derived neuronal cultures by reconstructing their effective connectivity, thus facilitating future efforts to generate predictive models for neurological disorders, drug testing and neuronal network modeling.

Summary

1. Introduction

• In-vitro cultures of primary neurons can self-organize into networks that generate spontaneous patterns of activity [10, 74, 58], in some cases resembling aspects of developing brain circuits [37, 24].
• These methods estimate the effective connectivity of a measured neuronal network by explicitly modeling the data generation process, i.e. only the connectivity of a simulated network model is inferred without any theoretical guarantee about its accuracy and its ability to correctly estimate the connectivity of the biological network [14, 76].
• Moreover, as the partial correlation matrix is symmetric, this method is not useful for detecting the causal direction of neuronal links.
• The authors show an experimental application of their approach to spontaneously generated in vitro networks of human iPSC-derived neurons cultured on MEAs providing an analysis and interpretation of the physiology not possible otherwise.

2. Materials and Methods

• The core of their methodology is graphically described in Figure 2 and includes three main phases: 1. statistical, correlationbased reconstruction of functional connectivity; 2. mathematically-rigorous super-selection of direct links via identification of peaks related to indirect and apparent links and 3.
• The algorithm constructs correlation triangles by considering all possible combinations of correlation delays for any possible triplet of neurons .
• The super-selection rule is formally presented later.
• Importantly, the algorithm removes correlations from the analysis, but does not remove inferred physical connections.

Temporal correlations

• (2) Using the Fast Fourier Transform (FFT) and the correlation theorem, computing correlations in Eq. (1) can be efficiently performed in O(S log(S)) with S the number of samples composing the signals sj .
• Their amplitude can be regarded as a measure of the level of correlation between the spikes in their registered firing activities.
• The higher the amplitude of a peak in Rjk(τ) (if any), the higher the probability that there is a statistical dependency between neuron j and neuron k.
• As explained earlier in the text, the existence of a functional coupling between two neurons does not necessarily imply that there is an effective connection between them.
• The location of each correlation peak with respect to the origin indicates the temporal delay τpeakjk between the activity of neuron j and neuron k .

Peak detection

• The authors implemented a peak detection algorithm applied to Rjk(τ) to identify all existing functional correlations between any pair of neurons (j, k) in the network and to discern the directional dependency between their spiking activities.
• As introduced above, the authors assume that a correlation peak in Rjk(τ) represents a potential connection between j and k and that τpeakjk is the signal propagation delay between them.
• The definition of “near-perfect match” is based on the choice for the parameter which represents the acceptable degree of temporal approximation when computing Eq. 5.
• The algorithm aims to construct njk × nkm × nmj correlations triangles given by all possible combinations of correlation peaks of any triplet of neurons’ pairs .
• Within this framework, each correlation triangle shares up to three edges with G, and the union E′ of all τhjk defines the direct graph G ′ = (V,E′).

Edge covering minimization

• Therefore, the challenge is to find a discrimination rule that allows us to select the correct τhjk in the correlation triangles that satisfy Eq. (4) (or Eq. (5)).
• The larger the variation of the delayed timing between the input and the output neurons, the larger the phase noise in the associated correlation peak which will have smaller amplitude and larger standard deviation than the correlation among neurons with input and output signals without phase noise and perfectly synchronized.
• To clarify this point, let’s consider the example reported in Figure 4.
• If the authors consider the triangles ACD or ABD as in Figure 4, in both cases, the delay τDA is the one detected and discarded because, as explained in Figure 3, the corresponding correlation peak has lower amplitude and larger width as it results from a signal that has propagated through three links.
• Within this picture, the authors can establish the second step in the super-selection algorithm as: Given a correlation triangle (τhjk, τ h′ km, τ h′′ mj), the delay associated with the smallest correlation peak am- plitude min(Ahjk, A h′ km, A h′′ mj) (8) estimates the false positive connection and is discarded from E′.

Connectivity matrix reconstruction

• There are three fundamental parameters that affect the inference performance of this method: T defines the time window (−T,+T ) over which to search for correlation peaks in the correlation function Rjk. σ is important because the location of the detected correlation peaks weakly depends on it.
• Large values increase the number of detected correlation peaks and computed correlation triangles with the potential shortcoming that true positive connections are filtered out.
• This process allowed us to individuate a range of values for T and σ where the inference performances of the algorithm reached a maximum plateau (Fig. 6D), and define the reconstruction rule of their connectivity method as follows (the numerical details are discussed in the next section of Materials and Methods).
• On the other hand, the maximum value for T can be chosen as several times (for example 5-6) the average delay among neurons in order to include more correlation triangles later used to refine the selection.

Experimental methodology

• To develop and validate their connectivity method the authors used spiking data generated via simulations of neuronal networks based on the Izhikevich model [29] (see Supplementary Material and Fig. S2).
• The original code was modified to guarantee high levels of activity in the network as well as bursting like behavior similar to that registered in their experiments .
• It is worth noting that their model is general and not restricted to sparse connectivity.
• For each network size, 20 different networks were randomly generated, simulated and then analyzed for connectivity reconstruction.

Performance measures

• The performances were assessed based on the indicators TP/Nc (net number of true positive (TP) connections, with Nc the number of connections in the simulated neuronal network), FP/Nc (net number of false positive (FP) connections) and ∆ = (TP − FP )/Nc (confidence indicator).
• ∆ is independent on the connectivity of the network being reconstructed and, because of its definition, it better defines the level of confidence in the detection of TPs by highlighting the method capabilities in rejecting or not the false positive connections.
• For the sake of comparison with the literature, the authors also used the more standard Accuracy measure ACC = (TP + FP )/Ntot [49, 19].

Numerical experiments

• This allowed us to isolate the activity of a few units per electrode, which resulted in the reconstruction of the activity of multiple detected neurons in the MEA well .
• In a few cases, the clustering approach did not perform efficiently and generated either too many or too few clusters which were detected by observing the increased number of outliers in the spike sorting output.
• For σ, observations on the recordings and the level of smoothing required for high-frequency noise removal led us to select σ ∈ (0.3, 0.8) ms.
• The authors calculated the average in-degree and out-degree of the vertexes of each directed graph corresponding to a reconstructed network by calculating all incoming and outgoing connections for each network vertex and then averaging over the total number of vertexes per network.

Numerical results

• To evaluate the performances of their method, the authors built a numerical model by designing and simulating an integrate-and-fire neuronal network mimicking the activity of in vitro cultures of neurons (see Supplementary Material).
• As a result, the same false positive connections observed for a given pair of neurons are very unlikely detected by many different points, and at high frequency they are filtered out.
• For two given points p1 and p2 the algorithm detected 15% and 20% FPs, respectively; however, most of the FPs detected by point p1 did not correspond to the ones detected by point p2.
• Since the authors made the hypothesis that each node in the network model has the same average connectivity (preserved in- and out-degree) which does not increase with the network’s size, they do not expect and they are not interested in seeing connectivity dependent variations.
• The authors also studied the integration and segregation properties of cultured human iPSC-derived neural networks (9C).

Discussion

• The authors demonstrated a new model-free based approach to infer effective - active, direct and causal - connections from in vitro neuronal networks recorded on MEAs.
• In terms of scalability, their algorithm matches the scalability properties of most of the state-of-the-art connectivity methods which, being based on pairwise statistical and correlation studies, scale with a practical computational complexity of the order of O(n2), where n is the number of neurons.
• More specifically, the authors have three main routines whose complexity should be assessed: the computation of pairwise correlations, the peak detection algorithm, and the detection of false positive connections.
• When applied to biological networks, their reconstruction method assumes that the spike sorting is accurate enough to separate the most prominent contributions in the activity of the network.
• High density MEAs that include tens of thousands of electrodes at cellular and subcellular resolution [73, 80, 38, 7], and complementary voltage and calcium imaging approaches [26, 41, 3, 27] that provide useful information about the localization and classification of recorded neurons, offer optimized acquisition settings and will greatly improve the resolution and accuracy of their approach.

4. Conclusions

• The authors have presented an innovative approach to map the effective connectivity of neural networks from multielectrode array data.
• Notably, their connectivity algorithm succeeds in detecting direct connections between neurons through a mathematically rigorous selection scheme that distinguishes between apparent or non-direct links and direct ones, therefore enabling inference of directed causal relationships between connected neurons.
• Furthermore, it will be broadly applicable to experimental techniques for neural activation and recording, increasing its utility for the analyses of spontaneous neural activity patterns, as well as neuronal responses to pharmacological perturbations and electrical and optogenetic stimulations [70, 26, 38, 41, 3].
• The authors thank Fabio L. Traversa for fruitful and constructive discussions on the theoretical framework.
• D.P. and A.B. designed and organized the experimental studies.

Figures

Figure 1: Definition of functional and effective connectivity. (A) Classification of causal (directional link), indirect (multi-neurons pathway) and apparent (functional coupling due to common input) connectivity. (B) Classification of most common connectivity inference methods in terms of causality and detection of direct links. On the x -axes, the graph shows a scale of causality which refers to the ability of a given connectivity method to infer or not the directionality of the functional connections between neurons. On the y-axes, the graph visually quantifies the capabilities of one approach to detect direct links between neurons by identifying and discarding multi-neurons connections and apparent ones. Indicators that are acausal and do not infer direct links can only report about functional connectivity (light yellow); indicators which contain information about direction of interaction and direct neuron-to-neuron communication are close to the inference of effective connectivity (orange color). Most common model-free techniques are indicated with black dots. Model-based methods are reported with blue triangles and their ability to infer effective connectivity is negatively weighted by the impossibility to test their performances. The super-selective correlation approach we propose is reported in black like the other model-free methods; a diamond signal is used to emphasize the fact that, although being model-free, it aims at inferring effective connections. The graph visually summarizes results of comparisons between connectivity methods from review papers [14, 76, 19] and do not contain precise quantitative information about the differences.

Figure 5: Real case scenario of correlation triangles and detection of false positive connections. (A) Correlation functions R1,4, R1,12, R4,12 computed for pairs of biological neurons (1, 4), (1, 12), (4, 12) (blue line) connected as in panel B. The selected time window was (-20 ms, +20 ms), but the graph only shows a zoomed view in (-10 ms, +10 ms) to better describe the selection rule if only one correlation was detected in time between the neurons. The σ of the Gaussian filter was 0.0005 s, resulting in smoothed correlation signals (orange dashed line). The detected correlation peaks (red stars) represent functional correlations between neurons and their combinations define the correlation triangles that our algorithm uses to detect false positive connections. The horizontal offset between the position of each peak and the origin represents the temporal delay τj,k associated with these links. (C) These delays satisfy the cyclic condition on τ , i.e., |τ1,4+τ4,12+τ12,1| = |τ1,4+τ4,12−τ1,12| < 3 ms therefore indicating that one of the three correlations in the correlation triangle corresponds to an apparent link. The algorithm detects which of the three by searching for the correlation peak with smallest amplitude. In the example, A1,4 has the smallest amplitude and the algorithm discards it from the analysis preventing generation of a false positive link.

Figure 4: Correlation triangles reconstruct all possible dependencies between neurons. On the left, a graph composed of four neurons A, B, C and D (polygon ABCD). Each edge is associated with a correlation peak and corresponding correlation delay τAB , τBC , τCD, and τDA. True connections are indicated by solid lines. The delay τDA corresponds to an indirect connection which the algorithm must identify and discard. This is indicated by a dotted line. On the right, the polygon ABCD is decomposed in four triangles ABC, ACD, ABD, and BCD (each color represents a different triplet). The participating correlation triangles reconstruct all possible dependencies among all neurons in the polygon ABCD. Each correlation link is associated with a correlation peak whose amplitude and width depends on the length of the indirect correlation path connecting the two correlated neurons (see Figure 3). In triangles ACD and ABD, the delay τDA results from a signal that has propagated through three links AB, BC and CD. Thus, τDA will be detected and discarded because having lower amplitude and wider peak than the other peaks.

Figure 6: Connectivity reconstruction. (A-C) Evaluation of the performance of the connectivity method under varying time window (−T,+T ), standard deviation σ, and for a fixed = 0.7 ms, for one simulated network of 10 neurons. (A) Net number of detected true positives TP/Nc. (B) Net number of false positives FP/Nc. (C) Ratio ∆ = (TP − FP )/Nc as chosen metric for evaluation. Nc is the known total number of connections in the simulated network. (D) Visual highlight on the behavior of ∆: a peak of performance is reached around 65%; low variability is proven for a wide range of T and σ (plateau indicated by the red plane). (E) Definition of a collection of K = 9 points p in the T -σ space for which the algorithm shows performances falling in correspondence of the plateau area. (F) Abstract representation of the statistical method for recognition of false positives connections and connectivity reconstruction. For each point a connectivity matrix Mp is computed, resulting in the computation of K = 9 different connectivity matrices for the same input network. Combination of the results enables computation of the frequency fjk for each connection following Eq. 9.

Figure 7: Method validation and performance evaluation. (A) Average (N=20) performances computed at different discrimination threshold d. The TPs remain roughly constant. On the other hand, the FPs decrease and fluctuate around very small percentage. The algorithm filters out the FPs that fluctuate at high frequency reaching best performances (85%) at d = 1. (B) Analysis of scaling properties. An average (N=20) number of TP/Nc, FP/Nc, and ∆ was computed for 20 randomly generated networks with 10, 20, and 50 neurons, respectively. Very good performances are maintained constant for increasing network size. (C) Accuracy ACC indicator computed for different network sizes as a function of the discrimination threshold d. The error bars stand for the standard deviation of a dataset of 20 different networks.

Figure 8: Connectivity reconstruction from neural recordings of developing hiPSC-derived neurons (A) Raster plots of the same MEA well at different time points during development: week 1, 2, 3, and 4 after plating. The four panels shows spiking signals from individual neurons (rows) obtained through spike sorting, PCA and k-means clustering of 300 sec recordings. The culture develops complex network features: from weakly active and randomly organized (individual spiking events), to very active and fully organized (network bursts). (B) Estimated effective connectivity of the developing culture whose activity is described in panel A. Each visual map consists of a 1.2 x 1.2 mm MEA plate (gray area), a 4-by-4 array of micro-electrodes (red circles) and the estimated directed connections (black arrows). The active neurons are represented by black dots; they are randomly distributed around their corresponding sensing electrode within a radius of 50 µm. (C) Directed graph relative to the culture at week 1 and week 4. The connectivity is equivalent to the one visualized in B but, for clarity and consistency with the main text, links between neurons are directed edges (arrows), active neurons are network nodes (blue: connected; yellow: independent). Given a MEA well and a specific time point, indexes refer to active neurons with recorded activity reported in (A). Note, neurons mapped at week 1 do not correspond to neurons mapped in the following weeks, although they are indicated with the same indexes.

Figure 3: The effect of indirect connectivity on correlation peaks. The edge AB, BC, CD, and DE are directed. The associated correlation peaks are high and narrow because there is a direct dependency between the activity of neuron A(B,C,D or E) and the one of neuron B(C,D or E). The links AC, AD and AE are indirect, with increasing order of correlation between the neurons because of more intermediary cells on the paths that link A to C, D and E. The longer the path, the weaker the dependency between the neurons which results in a lower and wider correlation peak.

Figure 9: Graph-theory based analysis of network connectivity in developing hiPSC-derived neuronal-network (A) Average number (N=20) of active neurons (black) and detected connections (red) in MEA wells recorded at week 1, 2, 3, and 4 after plating. The error bars represent the standard deviations for a dataset of 20 different MEA wells. (B) (Left panel) Statistics (N=20) on the neuronal in- (number of input connections) (blue) and out-degree (number of output connections) (yellow). The vertical red bars stand for the standard deviations of in- and out-degree data. (Right panel) Average maximum (continuous line) and minimum (dashed line) in- (top) and out- (bottom) degree. The number of in- and out- degree is equally distributed with equivalent raising behavior as a function of the cultures’ developmental stage. (C) Statistics on network hubs as a function of the culture’s age. In blue, average (N=20) maximum hub size. We used the in-degree centrality measures for hub detection and characterization. In yellow, percentage of neurons that function as network hubs relative to the total number of active neurons in the well (value averaged over 20 wells). More mature networks include few super-hubs (∼5% of the total number of neurons) with increase in size as demonstrated by the raising number of incoming connections at week 3 and 4. (D) Characterization of network segregation and integration properties. Values in the graph correspond to the mean Path Length (PL) and mean Clustering Coefficient (CCo) calculated for each well at week 1, 2 ,3, and 4 after plating and then averaged over 20 wells. The error bars correspond to the standard deviation of 20 different measures. The mean PL corresponds to the average shortest path length in the networks. Infinite (absent) connections between neurons were not considered in its calculation. The PL is very low for highly immature cultures where only sparse activity from individual neurons was mainly registered. It increases as soon as connections are formed (week 2) but no longer changes with the number of new connections (week 2, 3, 4). The CCo is low and decreases with the maturation of the network as indication of favored segregation versus integration with increasing number of independent but highly-integrated network units (hubs).

Figure 2: Connectivity reconstruction via detection of correlation triangles and classification of indirect and apparent links. 1. Given three neurons 1, 2 and 3, our algorithm searches, in time, for all correlations between them by computing the pairwise correlation functions R1,2, R2,3 and R3,1. In this representative example, the algorithm detects two correlation peaks for each pair of neurons and associates the corresponding delays of interactions (τ−11,2 , τ +1 1,2 , τ −1 2,3 , τ −2 2,3 , τ +1 3,1 , τ+23,1 ) which are defined by the location of the peaks with respect to the origin of the x axis. 8 possible combinations of interactions can occur in time between the three pairs of neurons. These correspond to the 8 correlation triangles shown on the right side of the panel. 2. Among the correlation triangles, the algorithm detects 2 critical cases (|τ−11,2 + τ −2 2,3 + τ +1 3,1 | < and |τ+11,2 + τ −1 2,3 + τ +2 3,1 | < ) and identifies the peaks relative to an indirect (multi-neurons path) and an apparent (common output) connection by searching for the ones with smaller amplitude. The smallest peaks (τ+13,1 and τ +2 3,1 ) are discarded from the analysis. 3. The correlation triangles are functional to the estimation of the direct connections in the network. Because no more correlation exists between neuron 1 and 3, the estimated effective connectivity includes only the direct links for (1, 2) and (3, 1): two connections with opposite directionality exist between neuron 1 and 2 because positive and negative correlation peaks are detected in R1,2 (τ −1 1,2 and τ +1 1,2 ); one link connects neuron 2 to neuron 3 as a result of the positive correlation peaks in R2,3 (τ +1 3,1 and τ +2 3,1 ).

Full text

Super-selective reconstruction of causal and direct connectivity with
application to in-vitro iPSC neuronal networks
Francesca Puppo
a,∗
, Deborah Pr´e
b,∗
, Anne Bang
b,∗∗
, Gabriel A. Silva
c,∗∗
a
BioCircuits Institute, Center for Engineered Natural Intelligence, University of California, San Diego, La Jolla, 92093 CA,
USA
b
Conrad Prebys Center for Chemical Genomics, Sanford Burnham Prebys Medical Discovery Institute, La Jolla, CA 92037,
USA
c
Department of Bioengineering, Department of Neurosciences, BioCircuits Institute, Center for Engineered Natural
Intelligence, University of California, San Diego, La Jolla, 92093 CA, USA
Abstract
Despite advancements in the development of cell-based in-vitro neuronal network models, the lack of ap-
propriate computational tools limits their analyses. Methods aimed at deciphering the effective connections
between neurons from extracellular spike recordings would increase utility of in-vitro local neural circuits,
especially for studies of human neural development and disease based on induced pluripotent stem cells
(hiPSC). Current techniques allow statistical inference of functional couplings in the network but are fun-
damentally unable to correctly identify indirect and apparent connections between neurons, generating re-
dundant maps with limited ability to model the causal dynamics of the network. In this paper, we describe
a novel mathematically rigorous, model-free method to map effective - direct and causal - connectivity of
neuronal networks from multi-electrode array data. The inference algorithm uses a combination of statistical
and deterministic indicators which, first, enables identification of all existing functional links in the network
and then, reconstructs the directed and causal connection diagram via a super-selective rule enabling highly
accurate classification of direct, indirect and apparent links. Our method can be generally applied to the
functional characterization of any in-vitro neuronal networks. Here, we show that, given its accuracy, it can
offer important insights into the functional development of in-vitro iPSC-derived neuronal cultures.
Keywords: neuronal network, effective connectivity, functional connectivity, apparent connectivity,
correlation, MEA, iPSC, development
1. Introduction
In-vitro cultures of primary neurons can self-organize into networks that generate spontaneous patterns
of activity [10, 74, 58], in some cases resembling aspects of developing brain circuits [37, 24]. The emer-
gent functional states exhibited by these neuronal ensembles have been the focus of attention for many
years [15, 83] as they can be used to investigate principles that govern their development and maintenance
∗
These two authors contributed equally
∗∗
Corresponding authors
Email addresses: email: abang@SBPdiscovery.org (Anne Bang), email: gsilva@ucsd.edu (Gabriel A. Silva)
Preprint submitted to Elsevier January 25, 2021
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprintthis version posted January 27, 2021. ; https://doi.org/10.1101/2020.04.28.067124doi: bioRxiv preprint

[35, 44] and to produce biological correlates for neural network modeling [11, 66]. The introduction of hu-
man induced-pluripotent stem cell (hiPSC) technologies [69, 82] opened the possibility to generate in-vitro
neuronal networks in typical [30, 42], as well as patient-specific genetic backgrounds [57, 5, 75, 78, 8, 39],
demonstrating the potential to reproduce key molecular and pathophysiological processes in highly con-
trolled, reduced, experimental models that enables the study of neurological disorders and the discovery and
testing of drugs, especially in the context of the individual patient [72, 16, 61].
One common approach to obtain information from in vitro neuronal networks is to record their activity
via multi-electrode array (MEA) or calcium fluorescence imaging and then use network activity features to
describe their physiology. One main limitation, however, is that these high-dimensional data, which report
about the information representation in the network, do not translate into a clear understanding of how this
representation was produced and how it emerged based on neuronal connectivity [14]. The synchronization
of spontaneous spike trains among different MEA sites or neurons, also referred to as network bursting, is an
example of observed neural behaviors widely reported in the literature. The generation of network bursting
in an in vitro neuronal culture is evidence that the neurons are synaptically connected. However, the extra-
cellular nature of the MEA recording does not provide information about how neurons are connected and
how signals propagate between them, such that computational analyses are necessary to reconstruct their
complex dynamic patterns and relate their emergence to the underlying wiring diagram [66]. However, this
kind of analysis presents several challenges as it requires not only identification of functional relationships
between cells, but also reconstruction of the dynamic causality (i.e., the knowledge of which neuron fires first
and affects another one) between directly linked neurons that are simultaneously involved in several differ-
ent signaling pathways. This defines the difference between functional and effective connectivity inference:
the first only reports about statistical dependencies between cells’ activities without giving any information
about specific causal and direct effects existing between two neurons [76]; the second attempts to capture a
network of effective - direct and causal - effects between neural elements [65].
Only model-based approaches [14, 34] have been proposed for inference of effective connectivity. Among
them, dynamic causal modeling (DCM) [18] and structural equation modeling [36] variants have shown best
performances. However, these methods estimate the effective connectivity of a measured neuronal network
by explicitly modeling the data generation process, i.e. only the connectivity of a simulated network model
is inferred without any theoretical guarantee about its accuracy and its ability to correctly estimate the
connectivity of the biological network [14, 76].
Because of this limitation, descriptive, model-free approaches are usually preferred as they are easy to
implement, rely on a limited number of assumptions that are directly related to the investigated neuronal
network, and can be more easily validated [14, 34]. A number of model-free methods proposed for recon-
structing the connectivity of in vitro neuronal networks [19] have been previously reviewed [48] and tested
2
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprintthis version posted January 27, 2021. ; https://doi.org/10.1101/2020.04.28.067124doi: bioRxiv preprint

[76]. However, because they rely on purely statistical indicators, they can only infer how neurons are func-
tionally coupled, but lack the ability to identify the network of effective interactions between neurons by
either missing the directionality or confounding indirect and apparent links from direct ones. Directional-
ity conveys the causality of signaling in the network, i.e. which neural element has causal influences over
another (Figure 1A(i)). However, causality does not imply a direct connection between two neurons. In
fact, a functional coupling between two neurons can be causal even though the two neurons are not directly
connected, and this may occur if there is a multi-neurons pathway between the two cells (indirect connec-
tion, Figure 1A(ii)), or if the connection detected between the two neurons is simply a mathematical artifact
resulting from the correlation of correlations generated by common inputs from other participating neurons
(apparent connection, Figure 1A(iii)) [17].
Methods such as correlation [53], coherence [25, 23, 12], mutual information [19, 22, 51], phase and
generalized synchronization [51, 2], and joint entropy [19] describe only statistical dependencies between
recorded neurons without carrying any information of causality or discriminating direct and indirect effects.
Techniques such as cross-correlation [19, 28], directed and partial directed coherence [1, 56], transfer en-
tropy [19, 28, 33] and Granger causality [21, 59] are examples of causal indicators as they provide inference
of directionality of dependence between time series based on time or phase shifts, or prediction measures.
However, because these operators rely only on pairwise statistical comparisons and treats pairs of neurons
independently, they show the same limitations when dealing with indirect connections and external inputs.
Only a few techniques can compete in the challenge of inferring the effective connectivity of a network.
Partial-correlation [19, 68], which takes into account all neurons in the network, showed best performance
in detecting direct associations between neurons and filtering out spurious ones [45]. The most significant
limitation of this solution is its high computational cost. Moreover, as the partial correlation matrix is
symmetric, this method is not useful for detecting the causal direction of neuronal links. It also does not
attempt to infer self-connections [14]. A combination of correlation and network deconvolution was used
by Magrans and Nowe [13] to infer a network of undirected connections with elimination of arbitrary path
lengths caused by indirect effects. However, this method also can not identify directions of connections and
the singular value decomposition of network deconvolution has an extremely high computational complexity
[45]. A convolutional neural network approach [54] showed the same limitations in computational complexity
and undetected self- and causal connections. Figure 1B graphically summarizes the inference capabilities of
the state-of-the-art connectivity methods as reported in [14, 76, 17, 2].
In this work, we propose a novel, mathematically rigorous method that uses a model-free approach (i.e.
does not depend on a set of underlying assumptions about the biology of participating cells) to decompose
the complex neural activity of a network into a set of numerically validated direct, causal dependencies
between the active component neurons that make up the network. First, the inference power of statistical
3
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprintthis version posted January 27, 2021. ; https://doi.org/10.1101/2020.04.28.067124doi: bioRxiv preprint

approaches (signal-, network- and information theory- based) allows mapping the functional connectivity
of the network. Then, we propose a mathematically rigorous selection scheme that distinguishes between
apparent or non-direct links and direct ones, therefore enabling inference of direct causal relationships be-
tween connected neurons which more realistically describe the effective connectivity of the network (Figure 2).
We evaluate the performances of the proposed method on synthetic datasets generated through simulation
of an integrate-and-fire neuronal network mimicking the activity of in vitro cultures of neurons, and demon-
strate important improvements, relative to the state-of-the-art connectivity methods, to network inference
accuracy due to a deterministic component of our method capable of identifying false positive connections.
We show an experimental application of our approach to spontaneously generated in vitro networks of
human iPSC-derived neurons cultured on MEAs providing an analysis and interpretation of the physiology
not possible otherwise. We describe the temporal evolution associated with the connectivity and dynamic
signaling of developing hiPSC-derived neuronal networks, including increasing synchronized activity and
the formation of small numbers of hyper-connected hub-like nodes, as similarly reported by others [30, 8].
These results further support the performance quality of our approach and provide an example of how this
connectivity method can be used to characterize network formation and dynamics, thus facilitating efforts
to generate predictive models for neurological disease, drug discovery and neural network modeling.
2. Materials and Methods
Theoretical framework for connectivity reconstruction
The central contribution of this manuscript is in providing an innovative, computationally efficient,
and easy-to-apply method for decomposing the collective firing properties stored in the electrophysiological
recordings from neuronal networks on MEAs into direct (one-to-one) and causal (directional) relationships
between all participating neurons. We propose a multi-phase approach that identifies and discards any
correlation link that does not directly relate to a direct interaction between two cells. The core of our
methodology is graphically described in Figure 2 and includes three main phases: 1. statistical, correlation-
based reconstruction of functional connectivity; 2. mathematically-rigorous super-selection of direct links
via identification of peaks related to indirect and apparent links and 3. reconstruction of directed causal
connectivity between neurons.
1. The functional connectivity (statistical dependencies) of a network is computed via pairwise correla-
tion studies. Functional interactions between neurons are represented by correlation peaks and their
delays τ . The algorithm constructs correlation triangles by considering all possible combinations of
correlation delays for any possible triplet of neurons (Figure 2.1). Importantly, correlation triangles do
not refer to any three-neuron physical connection, sometimes referred to as “neural triangles” in the
4
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprintthis version posted January 27, 2021. ; https://doi.org/10.1101/2020.04.28.067124doi: bioRxiv preprint

literature [63, 55]. Here we define a correlation triangle as a mathematical object that our algorithm
uses to classify functional interactions based on all possible triplets of correlation delays that can be
formed in the network. Therefore, correlation triangles exploit the entire signal history of neurons in
order to determine the correlation peaks.
2. Correlation peaks associated with indirect or apparent links in corresponding correlation triangles are
discarded from the analysis by means of a mathematical super-selection rule which deterministically
classifies the type of dependence between each triplets of neurons (Figure 2.2). The super-selection
rule is formally presented later. Here, it can be summarized as follows. If a correlation triangle is
made up of three correlation delays that are the combination of one another, one of the component
correlation delays is either representative of an indirect link (Figure 1A(ii)) or of an apparent link
(Figure 1A(iii)); therefore, this correlation delay does not refer to an effective connection and must be
discarded. When the algorithm finds a correlation triangle which satisfies this condition, it deepens
into the classification of the involved correlation delays and selects the correlation peak to remove
based on the peak’s amplitude. For example, in Figure 2.2, the algorithm identifies an indirect link
between 1 and 3 (a multi-neuron pathway), and an apparent link between 1 and 3 (correlation due to
a common output). The correlation peaks corresponding to these links in the correlation triangles are
discarded from the analysis. Importantly, the algorithm removes correlations from the analysis, but
does not remove inferred physical connections.
3. Only when all correlation peaks between two neurons are discarded, the algorithm recognizes that a
specific interaction is only apparent and deletes the corresponding connection. For example, in Fig-
ure 2.3, there is no existing connection between neurons 1 and 3. The estimated effective connectivity
includes only direct links for (1, 2) and (2, 3): two connections with opposite directionality exist be-
tween neuron 1 and 2 because positive and negative correlation peaks are detected in R
1,2
(τ
−1
1,2
and
τ
+1
1,2
); one link connects neuron 2 to neuron 3 as a result of the positive correlation peaks in R
2,3
(τ
+1
2,3
and τ
+2
2,3
).
The following sections describe the mathematical details of the developed technique. The connectivity
reconstruction algorithm and associated functions for performance evaluation were implemented in Matlab
and code is available online at https://github.com/fpuppo/ECRtools.git.
Reconstruction of functional connectivity
Temporal correlations
To identify the temporal correlations between the activity of all pairs of N recorded neurons j, k ∈
{1, ..., N} in the network, we computed the pairwise correlation function between the corresponding signals
s
j
and s
k
R
jk
(τ) =
Z
s
j
(t)s
k
(t + τ )dt (1)
5
(which was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
The copyright holder for this preprintthis version posted January 27, 2021. ; https://doi.org/10.1101/2020.04.28.067124doi: bioRxiv preprint

Frequently asked questions

Q1. What are the contributions in "Super-selective reconstruction of causal and direct connectivity with application to in-vitro ipsc neuronal networks" ?

In this paper, the authors describe a novel mathematically rigorous, model-free method to map effective direct and causal connectivity of neuronal networks from multi-electrode array data. Here, the authors show that, given its accuracy, it can offer important insights into the functional development of in-vitro iPSC-derived neuronal cultures. 

Q2. What have the authors stated for future works in "Super-selective reconstruction of causal and direct connectivity with application to in-vitro ipsc neuronal networks" ?

In addition, it has good scaling capabilities and can be further generalized to any kind of network, thus allowing to target different problems in intact neurons, synthetic models as well as in vitro and in vivo systems. Furthermore, it will be broadly applicable to experimental techniques for neural activation and recording, increasing its utility for the analyses of spontaneous neural activity patterns, as well as neuronal responses to pharmacological perturbations and electrical and optogenetic stimulations [ 70, 26, 38, 41, 3 ]. As novel electrophysiology technologies come online and are validated, the methods the authors presented here will be in an immediate position to take advantage of them, resulting in fundamental improvements in spatial resolution and reconstruction accuracy. Furthermore, their algorithm can be further extended, improved, and possibly integrated with already in-use techniques to overcome important limitations such as the detection of inhibitory connections and the inference of effective connectivity in the bursting regime.