Homophilic wiring principles underpin neuronal network topology in vitro
- MRC Cognition and Brain Sciences Unit, University of Cambridge, United Kingdom
- Department of Electrical and Electronic Engineering, Imperial College London, United Kingdom
- I-X, Imperial College London, United Kingdom
- Department of Physiology Development and Neuroscience, University of Cambridge, United Kingdom
- Department of Applied Mathematics and Theoretical Physics, University of Cambridge, United Kingdom
- Department of Biosystems Science and Engineering in Basel, ETH Zurich, Switzerland
- NRD, Roche Innovation Center Basel, F. Hoffmann-La Roche Ltd, Switzerland
- Department of Psychiatry, University of Cambridge, United Kingdom
- Division of Cognitive and Behavioral Neurology, Brigham & Women’s Hospital, United States
- Harvard Medical School, United States
Figures
Figure 1 with 2 supplements
High-density microelectrode array extracellular recordings of developing neuronal cultures.
(a) Scanning electron microscope (SEM) image of primary neurons plated on a high-density microelectrode array (HD-MEA; neurons are colored in red, electrodes in green; the electrode center-to-center distance is 17.5 μm). (b) Activity scan across the entire HD-MEA for an example sparse primary cortical (PC) neuron culture at DIV 14. Colors depict the average absolute amplitude of online-detected multi-unit activity per electrode (yellow colors indicate the location of potential neurons, i.e., high amplitude deflection, on the array; the array dimensions are 120×220 electrodes; electrode amplitude values are averaged over 1 min recordings). (c) Example extracellular trace recorded from one electrode (10 s), with high action potential (AP) spiking activity and bursts (middle panel); a single AP is depicted in the lower panel. (d) Electrical footprint (EF) of a single spike-sorted unit on the HD-MEA. The EF is the AP-triggered extracellular potential signature of a single unit on the HD-MEA, here depicted across 16 electrodes of a high-density recording block in light gray: 100 AP-triggered cutouts for this EF; in black: the average EF. The lower panel shows a spike autocorrelogram for this unit. (e) Tracking of individual EFs across development in vitro. The upper panel depicts the EFs inferred for the four recording time points (DIV 7, 10, 12, 14); the lower panel shows the average activity of the tracked unit (bin size: 100 s; gray to black colors correspond to the four recording time points; 30 min per recording day). (f) Spontaneous electrical activity of neuronal cultures, their average firing rate, increased significantly with development (n=6 sparse rodent PC networks). (g) A spike raster plot shows the spike-sorted activity for one culture (300 s-long zoom in on a network recording with 134 units); the lower panel depicts the binned activity of the same recording (bin size: 0.1 s). All neuronal networks showed a mixture of regular and more synchronized activity periods (network bursts). (h) In order to probe neuronal network development, we used the spike time tiling coefficient (STTC) to infer functional connectivity statistically. The four parameters required to calculate STTC connectivity graphs between Unit A and Unit B are PA, PB, TA, and TB. TA is calculated as the fraction of the total recording time ‘tiled’ by ±Δt of any spike from Unit A. PA is the proportion of spikes from Unit A, which lies within ±Δt of any spike(s) of Unit B (spikes in red). The spike trains for Unit A and B are depicted as black bars; ±Δt is depicted as blue lines for each spike. The significance of each STTC value was estimated against jittered surrogate spike train data. (i) The average STTC increased significantly with development (n=6 sparse rodent PC networks). (j) Example functional connectivity graph inferred from a DIV 14 PC network. Only the top 2% strongest connections are displayed; each dot represents the physical location of a putative neuron on the HD-MEA; the dot size corresponds to the nodal degree of the neuronal unit.
Figure 1—figure supplement 1
Tracking neuronal networks at cellular resolution on high-density microelectrode arrays.
Boxplots depicting the number of units exceeding a firing rate of 0.01 Hz in the sparse rodent network dataset (left panel; plating density 50,000 cells per well; n=6 cultures), and in the dense rodent network dataset (right panel; plating density 100,000 cells per well; n=12 cultures) over developmental time. Blue-colored boxplots/dots indicate the results obtained from the spike-sorting/postprocessing pipeline used in the main report (see Methods); the boxplots/dots in black show the alternative postprocessing pipeline (Kilosort 2.5 and Bombcell quality control; see Methods and Figure 7—figure supplement 2).
Figure 1—figure supplement 2
Global topological measures of sparse rodent cultures over time.
(a) Panel depicts the statistically inferred network density for the n=6 sparse rodent networks until days in vitro (DIV) 14. (b) Proportion of extant connection by distance and grouped by DIV 7 (pink), DIV 10 (light purple), DIV 12 (dark purple), and DIV 14 (dark blue). (c) The total degree for sparse rodent networks across development.
Figure 2 with 1 supplement
Probing wiring principles via generative network models.
(a) Example networks highlighting four common graph theoretical metrics. For each schematic, the node size corresponds to the respective graph statistic, and we provide the statistic within the node. The panel on the left shows a schematic for the degree (blue), which relates to the total number of edges that each node has. The left-middle panel shows the clustering coefficient (green), which is a measure of how nodes in a graph tend to cluster together; the right-middle panel shows the betweenness centrality (pink), which relates to the number of shortest paths in the network which pass through the node; the panel on the right shows a schematic for the total edge length (yellow), which relates to the total sum of the edge lengths among all its connections for each node. (b) Functional connectivity graphs, as inferred from high-density microelectrode array (HD-MEA) network recordings, were characterized by graph theoretical means. Each panel on the top row shows an example schematic network, with the node size corresponding to the respective graph statistic (degree, clustering, betweenness centrality, total edge length). The lower row shows histograms with the distribution of the metric across the network. (c) The generative network model works by simulating network development probabilistically according to the costs incurred in the network (here reflected by the Di,j term, which is modeled as the Euclidean distances) and the topological values. Both terms are parameterized, meaning that for each simulation, there is an altered tuning in how costs and values influence the probability with which nodes form connections iteratively. Parameters (η and γ) alter direction and magnitude to which the costs and values, respectively, influence wiring probabilities. For any combination of wiring parameters, we simulate a network developing over time until it reaches the number of connections, m, equal to the observed empirical network. (d) Comparing the simulated and empirically observed network, we fit the generative network model to achieve a simulation that is the least dissimilar to the observation. This is quantified by the energy equation and is shown by dark blue in the parameter space. The energy equation is given by the maximum Kolmogorov-Smirnov (KS) statistic across degree, clustering, betweenness centrality, and edge length. Each parameter combination is plotted on an energy landscape, which demarcates the model fits across the two-dimensional parameter space. Lower energy values correspond to better fits.
Figure 2—figure supplement 1
Distribution of inter-neuronal Euclidean distances across rodent and human datasets.
(a) Density plot of inter-neuronal Euclidean distances for each culture. Note that each matrix was used as the Di,j term for the generative model constructed for each time point of that culture. For all distributions, we plot the total set of distances between neurons, irrespective of the inferred functional network. Panel a depicts Euclidean distances for the sparse rodent networks (50,000 cells per well); the top panel shows the overall distributions; the bottom panel shows the individual distance matrices. (b) Euclidean distances for the dense rodent networks (100,000 cells per well). Panels (c-e) show Euclidean distance distribution across the iPSC-derived human neuron lines at DIV 28 (c, glutamatergic; d motor and e dopaminergic neurons). (f) Euclidean distances for the human cerebral organoid recordings.
Figure 3 with 5 supplements
Homophilic principles underpin rodent primary network development and explain variance in their maturational trajectories.
(a) The top row shows a representative sparse rodent primary cortical (PC) network developing over the first two weeks in vitro on the high-density microelectrode array (HD-MEA) (DIVs 7, 10, 12, and 14; sparse neuron plating density; gray nodes correspond to single neurons; the size of each neuron corresponds to its nodal degree). Below, we show the energy of the top performing (n=1) simulated networks, given across the four main categories of studied wiring rules (with all models plotted, grouped by their rule). Boxplots represent the median and interquartile range (IQR); outliers are demarcated as small black crosses and are those which exceed 1.5 x the IQR away from the top or bottom of the box. (b) Exemplary immunohistochemical staining of a rodent PC culture (control experiment; NeuN staining in blue, βIII-tubulin in red, and Synaptophysin in green). (c) All energy landscapes of DIV 14 networks, plotted on top of each other, under the matching generative network model (a homophily model). (d) Cumulative density functions (CDFs) of the matching generative model, showing that simulated and observed statistics overlap very well. CDFs are shown across the four network statistics included in the energy equation of the top = 99 simulations for the matching model (best performing generative model) compared to an observed network (black line). Subsequent prank values were computed using a Monte Carlo sampling procedure. (e and f) Simulated network trajectories were examined to determine if the later generative model’s trajectory (DIV 14) recapitulated earlier observed statistics. Simulated network statistics were derived by computing the modularity Q statistics (e) and global efficiency (f), at each time-point, within the generative model that was best fit to the DIV 14 networks. Results were scaled to the age of the culture so that observations and simulations could be compared directly, where 50% of the simulated time corresponds to DIV 7 and 100% of simulated time corresponds to DIV 14. Note that simulated time here refers to the percentage completion of the simulation, rather than the number of connections added, so that cultures could be aligned based on time. The quoted r and p values correspond to the Pearson’s correlation between observed and simulated statistics shown.
Figure 3—figure supplement 1
Generative network modeling results for the sparse rodent cultures.
Generative model fits (energy), including all 13 wiring models, for the sparse rodent networks (50,000 cells per well plating density; n=6 cultures) across development. Each boxplot presents the median and interquartile range (IQR). Outliers are demarcated as small black crosses and are those which exceed 1.5 x the interquartile range away from the top or bottom of the box. Generative model performance over time according to the energy equation. In each box, the energy of the top a n=1, b n=10, and c n=50 performing simulations, respectively.
Figure 3—figure supplement 2
Generate network model fits as a function of binarization threshold.
Comparison of (a) network density and (b) generative network model fits for the sparse rodent networks (50,000 cell per well plating density; n=6 cultures) at DIV 14 indicates that the observed ordering in energy values is robust against different significance thresholds. As expected, the network density decreased with more conservative thresholds (smaller alpha values). For all alpha values <0.05, we can reproduce the same ordering of network model fits as reported in the main manuscript. Boxplot colors and ordering of generative models are the same as depicted in Figure 3—figure supplement 1. Each boxplot presents the median and interquartile range (IQR). Outliers are demarcated as small black crosses and are those which exceed 1.5 x the interquartile range away from the top or bottom of the box.
Figure 3—figure supplement 3
Generative model outcomes as a function of empirical network data.
(a) For the sparse rodent networks (50,000 cells per well; n=6 cultures × 4 time-points) across development (DIV 7, DIV 10, DIV 12, DIV 14) we computed the Pearson correlation between the empirical activity/network measures (x-axis) and the energy (left box), η (middle box), γ (right box) across generative models (key provided on the right). Only significant correlations are provided, with non-significant (p>0.05) values grayed out. The circle, triangle and star shapes highlight what plots are plotted subsequently in the remaining panels and pertain to the matching model (generative model 2) and the best performing non-homophily model - the degree average model (generative model 3).(b) Under the matching generative model, the relationship between model energy and the mean spike time tiling coefficient (STTC) (left), network density (middle) and modularity (right) is non-significant. (c) Under the degree average generative model (the best non-homophily) model, the relationship between model energy and the mean STTC (left), network density (middle) and modularity (right) is positive. This suggests that as networks develop and form more correlated activity, become denser and less modular, this model performs worse. (d) Under the matching generative model, the relationship between the optimized η and the mean STTC (left), network density (middle) and modularity (right) is non-significant. (e) Under the matching generative model, the relationship between the optimized γ and the mean STTC (left), network density (middle) and modularity (right) is non-significant.
Figure 3—figure supplement 4
Alternative generative model procedures.
Generative model fits (energy), including all 13 wiring models, for sparse primary rodent networks (n=6 cultures) at DIV 14. (a) Generative models fit only using the wiring term (Ki,j) when forming networks. (b) Random null models. Completely randomized networks (size- and density-matched, left) and networks which have been partially randomized via the randmio_und() function, using the Maslov-Sneppen algorithm (Maslov and Sneppen, 2002), which rewires each edge approximately five times. In the case of completely random networks (left), there is no difference in relative performance across wiring rules, because the optimal parameter fit is achieved when both wiring parameters are set to zero, as this generates a totally random topology regardless of the rule. (c) Inferred transfer entropy (TE) networks. As the TE is a directed graph, we show the results for the symmetrized network (left); out-connections (middle) and in-connections (right). For each graph, the boxplot presents the median and interquartile range (IQR). Outliers are demarcated as small black crosses and are those which exceed 1.5 x the interquartile range away from the top or bottom of the box. Generative model performance over time according to the energy equation. In each box, the energy of top n=1 performing simulation is shown. (d) The relationship between the energy acquired from the same networks derived from the spike time tiling coefficient (STTC) (x-axis) and the TE algorithms (y-axis) is consistently positive; with TE derived from symmetrical TE connections, i.e., keeping only reciprocal edges (left; r=0.849, p=8.97 × 10-22), derived from all out TE connections (middle; r=0.867, p=9.55 × 10-25), derived from all in TE connections (right; r=0.901, p=2.91 × 10-29).
Figure 3—figure supplement 5
Generative model fits recapitulate observed network statistics.
Cumulative density functions (CDFs) of the top = 99 simulations (of the total 20,000 simulations) for the best performing generative models in each class (a, spatial. b, matching. c, clustering average. d, degree average) across the four statistics included in the energy equation (top four panels) and two additional metrics - the local efficiency (gray) and participation coefficient (yellow) - not included in the energy equation (bottom two panels). For visualization, we show only the solutions for a single sparse rodent culture. p values were computed using the Monte-Carlo bootstrapping procedure (see Methods; Generative network models).
Figure 4 with 2 supplements
Homophilic generative mechanisms capture the local topology in developing neuronal networks.
(a) Schematic illustrating the inference of topological fingerprints (TFs). A TF measure is computed as the Pearson’s correlation matrix among several topological statistics at the nodal level (degree, clustering, betweenness, edge length, efficiency, and matching). Each node in a corresponds to a single neuronal unit/neuron. The right panel shows the negative correlation between clustering and betweenness, which captures an aspect of the topological structure of the network shown on the left. The color bar is clipped to +/-1 (blue/yellow) for clarity. (b) The TFdissimilarity measures the extent to which GNM simulations capture the topological structure of the experimentally inferred networks. Homophily generative models, and spatial models show the lowest TFdissimilarity, suggesting that both can reconstruct local connectivity patterns in vitro (n=6, sparse PC networks; boxplots present the median and IQR; outliers are demarcated as small gray crosses, and are those which exceed 1.5 x the interquartile range, IQR). (c) Averaged TF matrix for the empirically observed data (on the left; DIV 14), versus the GNM results obtained from models with the best fits, i.e., the lowest energy values obtained for each model class: the matching rule (top left panel), the clustering-average rule (bottom left panel), the degree-average rule (top right panel) and the spatial mode (bottom right panel). Each node-wise measure used within the correlation matrix is plotted (left). As matching is an edge-wise measure, the matching calculation presented on row six was derived as a node-wise average. (d) This plot depicts the relationship between energy and TFdissimilarity values for each sparse PC network broken down by generative rule class (spatial in yellow, homophily in red, clustering in green, degree in blue). Each dot indicates the value of the two model fit functions for a single simulated network.
Figure 4—figure supplement 1
Homophilic generative mechanisms account for local relationships in developing dense rodent neuronal networks.
(a) Homophily generative models produce the lowest topological fingerprint (TF) dissimilarity for dense rodent networks (DIV 14 and 28), suggesting that it can reconstruct the local connectivity patterns of in vitro neuronal networks. In total, there are n=12 data points (one per culture) shown in each of the 13 boxplots. Boxplots present the median and interquartile range (IQR). Outliers are demarcated as small gray crosses and are those which exceed 1.5 x the interquartile range away from the top or bottom of the box. (b) A visualization of the averaged topological organization matrix for the observed (left), matching (middle left), degree-average (middle), clustering-average (middle right), and spatial (right) models at DIV 28.
Figure 4—figure supplement 2
Relationship between model energy and topological fingerprint dissimilarity in dense rodent neuronal cultures.
Using dense primary cortical (PC) rodent networks (100,000 cells per well), we show n=156 data points (n=12 cultures x n=13 generative model simulations) corresponding to the top n=1 performing simulation’s energy and its topological fingerprint dissimilarity performances at DIV 14.
Figure 5 with 2 supplements
Effect of plating density on network topology and generative network principles.
(a) Correlation matrix for global network statistics calculated across sparse (lower triangle) and dense (upper triangle) primary cortical (PC) cultures shows a highly similar covariance structure; colors indicate the Pearson’s correlation coefficient. Elements (5,9) and (9,5) of the correlation matrix are further highlighted in panel. (b) This plot shows the correlation between efficiency and small-worldness, for sparse and dense networks (r=−0.822, p=2.73 × 10–14); each point in the scatter plot corresponds to a single network. (c) Sparse (50,000 cells per well) and dense (100,000 cells per well) PC networks can be best simulated by the homophily generative models (DIV 14). Note that the leftmost boxplot is the same as that given in Figure 3a. (d) The energy landscapes for the matching generative network model landscapes inferred for the sparse and densely plated PC networks (DIV 14; see also Figure 3c).
Figure 5—figure supplement 1
Comparison of global network statistics across sparse and dense rodent networks at DIV 14.
(a) Comparisons were computed for network density, total local efficiency, edge length, matching, betweenness, modularity, strength, total degree, clustering and small-worldness. For each comparison, we provide the p-value computed from a Mann-Whitney U test. Sparse rodent networks were plated at 50,000 cells per well, dense networks were plated at 100,000 cells per well. Note here that the edge length refers to the average length of existing functional connections in the network. This measure contrasts with Figure 2—figure supplement 1 which shows the inter-neuronal Euclidean distances between all neurons, irrespective of the inferred functional network. (b) Direct comparisons of the generative model fits (energy), including all 13 wiring models, for the sparse (left) versus dense-plated (right) rodent networks at DIV 14. Boxplot colors and ordering of generative models are the same as depicted in Figure 3—figure supplement 1. The boxplot presents the median and interquartile range (IQR). Outliers are demarcated as small black crosses and are those which exceed 1.5 x the interquartile range away from the top or bottom of the box. Generative model performance over time according to the energy equation. In each box, the energy of top n=1 performing simulation is shown.
Figure 5—figure supplement 2
Generative network modeling results for the dense rodent cultures.
Generative model fits (energy), including all 13 wiring models, for the dense rodent networks (100,000 cells per well; n=12 cultures) across development showing the top-performing parameter combination. Boxplots present the median and interquartile range (IQR). Outliers are demarcated as small gray crosses and are those which exceed 1.5 x the interquartile range away from the top or bottom of the box.
Figure 6 with 4 supplements
Chronic block of GABAA receptors alters the wiring parameters of the homophily generative model.
(a) Representative neuronal population activity of a control/untreated (red, top) and gabazine-treated networks (bottom, blue) at DIV 14. In the gabazine condition, network burst activity is more strongly synchronized compared to the control networks. (b) Panel depicts that the number of putative inhibitory connections increased significantly following the washout of gabazine at DIV 14. Putative inhibitory connections were inferred by a cross-correlogram approach and using an algorithm to infer spike transmission probability, respectively, spike suppression (English et al., 2017). (c) The energy of the top-performing (n=1) simulated networks, given across the four main categories of studied wiring rules. Boxplots represent the median and interquartile range (IQR); outliers are demarcated as small black crosses and are those which exceed 1.5 x the IQR away from the top or bottom of the box. Energy values and wiring parameters for each individual rule of the top performing (n=1) simulated networks across control, gabazine-treated, and washout cultures are shown in Figure 6—figure supplement 4c, respectively. (d) The average probability score kernel density distributions (Pij) for the control (brown), gabazine (cyan), and washout conditions (lilac) recorded in sparse primary culture networks across simulated network development. The average gabazine probability distribution is shifted to the right and flatter compared to the controls and the washout condition. This indicates that wiring becomes more random across all time points. In Figure 6—figure supplement 4d we show all distributions used to construct these averages. Asterisk indicates p<0.05.
Figure 6—figure supplement 1
Chronic gabazine application to block GABAA receptor activity led to changes in spiking patterns.
(a) Representative spike train raster plots from a control (upper panel), gabazine-treated culture (middle) and washed out gabazine-treated culture; the panel below shows the representative population activity vectors (activity is aggregated 1 s bins). Panels on the right of the temporal raster plots show example close-ups of 1 s time windows taken from time bins comprising either median network activity (left, indicated by cyan triangle) or peak network activity (right, indicated by yellow triangle). (b) A logarithmic histogram of interspike intervals (ISIs; upper panel) showing that gabazine relatively increases the proportion of very short ISIs and long ISIs, giving a bi-modal distribution, equivalent to periods of relative quiescence and periods of fast spiking (bursts). The spike time tiling coefficient (STTC) distribution (probability plotted on a log10 scale for visualization) shows gabazine-treated cultures had a far greater number of strongly correlated units compared to controls in which there were far fewer strong functional connections. (c) Cultures treated with gabazine showed decreased mean firing rate (Hz) of cellular activity compared to controls, but increased and more regular spiking and bursting activity (reduced coefficient of variation in inter-spike and inter-burst intervals). Gabazine treatment also led to much stronger edge weights (STTC values) and less positively skewed edge weight distributions (less positive skewness of STTC distribution)—this skewness in controls reflects the presence of high strength (hub) nodes. Asterisks, *, ** and *** indicated p<0.05, 0.01, and 0.001, respectively.
Figure 6—figure supplement 2
Chronic gabazine application to block GABAA receptor activity led to changes in functional network topology.
(a) Representative examples of topological network plots from a sparse DIV 14 rodent culture (upper panel) and a culture treated with gabazine from DIV 1 (lower below). Each node (colored circles) represents the spatial location of a putative neuronal unit detected through spike sorting. Each edge (gray lines) represents a significant correlation in spiking activity between a pair of nodes. Functional connections were inferred using the spike timing tiling coefficient (STTC). The size of each node is proportional to the number of significant functional connections; its color reflects its spike rate. (b) Each dot represents a single neuronal unit. For both control cultures (upper; brown) and gabazine-treated cultures (lower; blue), there was no clear linear relationship between spike rate and mean STTC (average weight of all significant edges of a given node). Each group had similar firing rate distributions, though very different average functional connection weight distributions with gabazine-treated units being hyperconnected compared to controls. (c) A range of global topological metrics inferred from functional connectivity graphs of control (C) and gabazine (G) treated cultures. For each comparison, we provide the p-value computed from a Mann-Whitney U test. Asterisks, *, ** and *** indicated p<0.05, 0.01 and 0.001, respectively.
Figure 6—figure supplement 3
Inferring inhibitory connections using spike-transmission probabilities.
(a) An example deconvolved cross-correlation histogram (dCCH) to estimate the spike-transmission probability (STP) between two neurons in a culture treated with gabazine, and after washout. The median spike count prediction across time lag values is indicated by the orange line. (b) The same unit pair as in panel a, however, plotted with spike count relative to the prediction (dCCH after washout of gabazine). Red bars indicate where the spike count was lower than predicted in neuron j between 0–5ms after neuron i fired, thus indicating a putative inhibitory connection.
Figure 6—figure supplement 4
Comparison of generative model outcomes between control, gabazine-treated and washout conditions.
(a) Direct comparisons of the generative model fits (energy), including all 13 wiring models, for control (left) versus gabazine (middle) and washout (right) condition. The boxplot presents the median and interquartile range (IQR). Outliers are demarcated as small black crosses and are those which exceed 1.5 x the interquartile range away from the top or bottom of the box. In each box, the energy of top n=1 performing simulation is shown. (b) The energy landscapes of the matching generative rule for control (left) versus gabazine (middle), and after washout (right). (c) Wiring parameters of control (left; colored in red), gabazine (middle; colored in blue) and after washout (right; colored in purple) derived from the top n=1 performing matching simulation in terms of the wiring equation are shown. (d) The probability distributions for the control (top), gabazine (middle) and washout (bottom) condition were computed and scaled at 1% increments throughout the developmental course of the simulations. These distributions were averaged over time to form the mean probability distributions shown in Figure 6d.
Figure 7 with 3 supplements
Generative network modeling across human monolayer neuron cultures and cerebral organoids.
(a) Immunohistochemical staining of a control dopaminergic neuron (DN) network (top panel: expression of GFAP in red, MAP2 in blue, TH in green, and DAPI in gray). (b) Immunohistochemical staining of a single human cerebral organoid (hCO) slice (Tau in green, NeuN in purple and DAPI in blue; bottom panel: Tau in green, GFAP in purple, and DAPI in blue). (c) Clustering of human and rodent neuronal cultures based on a t-distributed stochastic neighborhood embedding (tSNE). Each dot corresponds to one culture; similarity was estimated by correlating the spike train autocorrelograms across all datasets and by then applying the tSNE. (d) Representative population activity plots for networks of glutamatergic (left), motor (left middle), dopaminergic (right middle) iPSC-derived neuronal cultures and a hCO slice (right). (e) Global topological measures of human iPSC-derived neuronal networks at DIV 28: the mean spike time tiling coefficient (STTC) (top left), network density (top right), small-worldness (bottom left) and proportion of extant connections by distance (bottom right). (f) STTC functional connectivity measures and small-worldness of hCO slices including: average STTC (top left), network density (top right), small-worldness (bottom left), and connectivity as a function of distance (bottom right). (g) Functional connectivity graph inferred from a hCO slice. (h) Fits across generative model rules in human iPSC-derived neuronal networks (DIV 28) and hCOs (right). The energy of the top n=1 performing simulated networks are shown. Boxplots represent the median and interquartile range (IQR). Outliers are demarcated as small black crosses and are those which exceed 1.5 x in the IQR away from the top or bottom of the box. The asterisk reflects homophilic rules, which demonstrate statistically lower energy (p<0.05) than other rules. All statistics are provided in Supplementary file 1h.
Figure 7—figure supplement 1
Spiking dynamics in human monolayer neuron and cerebral organoid cultures compared to rodent cultures.
(a) Temporal raster plots indicate the firing rates (spikes/s) of each cell line over consecutive 1 s time bins through 5 min of recording. Of the iPSC lines shown, motor neuron patterns most closely resemble the example rodent networks showing similar population bursts but lacking the variability across neurons, such as the subset of fast-spiking neurons (upper rows of rodent culture). 28-day-old dense rodent primary cortical (PC) networks (100,000 cells per well) are included here for comparison, as this is the latest time point recorded for the rodent PC cultures, and ostensibly the peak maturity in terms of network activity. Both glutamatergic and dopaminergic iPSC cultures showed highly regular network bursting activity, though at differing frequencies. Human cerebral organoid networks had fewer units but had periods of much higher firing rates. For each culture type, autocorrelograms of spiking activity are plotted to the right of each temporal raster plot. Autocorrelograms were constructed by first cumulatively concatenating spike times across all units in a recording into a single vector of spike times. The count of synchronous spikes across the network (in 1ms bins) was quantified across a range of time lags (zero lag was excluded). This indicates the degree of coactivation at various latencies. In all cell culture types, there was much higher spiking synchrony at shorter latencies indicated by peaks close to 0ms lag. However, rodent primary cultures (PCs) and human motor neurons (MNs) showed sharp declines in autocorrelation at longer latencies, indicating a strong preference for short latency, high frequency coactivation. Human glutamatergic (GNs) and dopaminergic (DNs) cultures had much greater autocorrelation at lag values further away from zero. This indicates coactivation at longer latencies and oscillatory network activity at slower frequencies. Human cerebral organoids (hCOs) were much more variable. (b) The autocorrelation vectors (autocorrelation at each lag value) across all cultures were correlated and a high level of concordance within culture sources in terms of latency.
Figure 7—figure supplement 2
Replication of the main results using an alternative spike-sorting and postprocessing pipeline.
(a) Panel depicts the generative model energy values for the sparse (1000 cells/mm2) rodent cultures across DIV 7–14 (n=6 cultures; functional connectivity inferred through spike time tiling coefficient (STTC)) using an alternative spike sorting and postprocessing pipeline (see Methods). Panel b shows the results for the dense (2000 cells/mm2) cultures, and panel c shows the results for the human neuronal networks (GNs: glutamatergic neurons; MNs: motor neurons; DNs: dopaminergic neurons); the alpha value to threshold the STTC matrices was 0.001. Panel d depicts the results of the TE analysis (alpha value to threshold TE matrices: 0.01). Since the used generative network modeling algorithms had been developed for symmetric binary connectivity matrices, we show results across five variants of differently symmetrized TE matrices (left: keeping only the reciprocal connections; middle left: keeping all connections; middle/middle right: symmetrizing either the in or out connections; right: proportional thresholding of the symmetrized weighted TE matrices and keeping only the strongest 5% connections). Boxplot colors as defined in Figure 3—figure supplement 1. Boxplots present the median and IQR. Outliers are demarcated as small gray crosses and are those which exceed 1.5 x the interquartile range away from the top or bottom of the box.
Figure 7—figure supplement 3
Replication of the main results on a validation dataset.
Generative network modeling results calculated on an independent dataset, comprising n=12 rodent networks at the sparse plating density (50,000 cells per well or 1000 cell/mm2). The recordings were performed at DIV 14. Results are based on binarized spike time tiling coefficient (STTC) functional connectivity matrices (alpha value to threshold STTC matrices: 0.001). The results of this validation dataset strongly resemble the original results reported for the sparse rodent networks in the main manuscript. Boxplot colors as defined in Figure 3—figure supplement 1. Boxplots presents the median and interquartile range (IQR). Outliers are demarcated as small gray crosses and are those which exceed 1.5 x the interquartile range away from the top or bottom of the box.
Tables
Table 1
Non-default parameter values used for the continuous-time TE inference.
| Parameter | Description | Value |
|---|---|---|
| lX | Number of interspike intervals in target history embeddings | 4 |
| lY | Number of interspike intervals in source history embeddings | 2 |
| kglobal | Number of nearest neighbours to find in the initial search | 10 |
| kperm | Number of nearest neighbours to consider during surrogate generation | 4 |
| Nsurrogates | Number of surrogates to generate for each node pair | 100 |
| sampling_method | Method with which to place the random sampling points. | "jittered_target" |
| jittered_sampling_noise | Width of the uniform jitter added to the target spike times used in resampling when sampling_method is set to "jittered_target". | 200.0 |
Additional files
-
Supplementary file 1
Supplementary tables.
(a) Overview of the datasets used in the study. (b) A list of all the value Ki,j terms that were included in the generative modeling, as given in the wiring equation. (c) A summary of all optimized parameters and energy values, for each dataset, across generative rules. (d) Statistical comparisons of the sparse rodent culture energy values across generative rules. (e) Statistical comparisons of the sparse rodent culture topological fingerprint dissimilarity across generative rules. (f) Statistical comparisons of the dense rodent culture topological fingerprint dissimilarity across generative rules. (g) Statistical comparisons of the dense rodent culture energy values across generative rules. (h) Statistical comparisons of human monolayer neuron and cerebral organoid culture energy values across generative rules. (i) Overview of used antibodies. (j) Confidence intervals of model energies, generated by subsampling from the dense rodent networks.
- https://cdn.elifesciences.org/articles/85300/elife-85300-supp1-v2.docx
-
MDAR checklist
- https://cdn.elifesciences.org/articles/85300/elife-85300-mdarchecklist1-v2.pdf
Download links
A two-part list of links to download the article, or parts of the article, in various formats.
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Homophilic wiring principles underpin neuronal network topology in vitro
eLife 14:e85300.
https://doi.org/10.7554/eLife.85300
