Reconstructing synaptic input activity during spontaneous network spiking

Schematic of main concept. Long-term whole-network extracellular recordings of neuronal spiking are performed using a high-density microelectrode array (HD-MEA) platform with 26’400 electrodes and 1024 channels for simultaneous readout. Next, the incoming monosynaptic connections to individual postsynaptic cells are identified using a novel regression approach applied to simultaneous HD-MEA and whole-cell patch-clamp recordings. Finally, spike trains (si) of all identified excitatory (PreE) or inhibitory (PreI) presynaptic cells are convolved with their respective postsynaptic-current estimate (wi) to obtain, for several hours, parallel reconstructions of the excitatory and inhibitory synaptic input activities. Crucially, postsynaptic spike times are also available during the reconstruction period. The experimental implementation of this approach is shown in Fig. 2. Main analytical steps are introduced in Fig. 3 and 4.

Experimental pipeline

(A) Left: Sensing area of a HD-MEA chip and magnified region with cultured primary neurons. Right: whole-array spike rate and amplitude maps. Electrodes displaying spiking activity were selected for subsequent long-term recording steps (grouped in 2-3 subsets of 1024 electrodes each).

(B) Left: example electrode traces of a long-term recording. Right: spike raster plot of one entire electrode subset. The total recording time for each electrode was at least 3 h.

(C) Example spike-triggered average extracellular signatures (or ‘footprints’) of individual units. A black asterisk marks the footprint magnified in the inset. Typically 100–200 units were identified per subset of 1024 electrodes.

(D) Left: fluorescence images of a patched cell. Middle: example traces from a simultaneous patch-clamp (green) and HD-MEA (blue) recording. Right: spike-triggered averaging of the extracellular signals, based on spike times detected via the patch-clamp electrode, revealed the HD-MEA footprint of the patched cell (blue traces). The footprint of the patched cell was matched to a unit footprint from the preceding spike-sorted long-term recording (orange).

(E) Electrode spike raster plot (top) and voltage-clamp trace (green; bottom) of a simultaneous HD-MEA and whole-cell patch-clamp recording. This paired recording was used to identify incoming monosynaptic connections (see Fig. 3).

In all panels, HD-MEA electrode traces were band-pass filtered at 0.3–9.5 kHz.

Connectivity inference and EPSC estimation based on paired HD-MEA and patch-clamp recordings

(A) Results from one example paired recording (e.g., Fig. 2E). We developed a linear regression procedure that estimates the EPSC that was evoked in the patched cell by each of the extracellularly recorded neurons. Middle: EPSCs of putative connections (amplitude > 10 × standard deviation of pre-spike baseline & amplitude > 1 pA). Black arrowhead indicates presynaptic spike time. Right: EPSC estimates of putatively unconnected units. Left: approximate spatial distribution of the presynaptic and postsynaptic cells. Locations of unconnected units are indicated by grey dots.

(B) For the same paired recording as in (A), example recording period with a raster plot of presynaptic unit spiking (top), the VC patch-clamp recording (black), and two current-trace reconstructions using unit spike times and EPSC estimates of either all units (magenta) or of the putatively connected units only (green).

(C) Variance of the measured patch-clamp recording explained by the respective current-trace reconstruction with all units (14 patched cells, with 142 identified connections, from 5 preparations). Different frequency contributions were compared (using fourth order Butterworth filters). Box plot indicates median and interquartile ranges, and whiskers the minimum/maximum values (excluding outliers).

(D) Validation of the regression-approach by simulation of ground-truth synaptic inputs. Left: comparison of two example EPSC estimates with the corresponding ground-truth EPSC. Right: mean errors between ground-truth and EPSC estimate (top) and F1 score (bottom) for different simulation parameters (n = 14 simulations each; see Methods for details).

See also Fig. S1 for evidence that the variation in the number of identified inputs is of biological origin.

Connection-type classification based on network-wide spike transmission or suppression

(A) Top: auto- (red, excitatory cell; blue, inhibitory cell) and cross-correlograms (CCG; black) for four example units from a long-term HD-MEA recording (> 3 h). The schematic indicates putative connections. Bottom: illustration of spike-transmission probability (STP) extraction. The dashed blue line (left) indicates the CCG baseline.

(B) Example matrix of STP values for all pairwise combinations of the units identified in a long-term HD-MEA recording. The rows labeled Ia/b and Ea/b correspond to the units shown in (A). The mean STP (across the sub-network) for each reference unit is shown on the right. Units are classified as being putatively excitatory or inhibitory based on the sign of their mean STP value. The color scale values in brackets apply to the column of mean STP values. Note the relatively large STP values associated with E-I compared to E-E unit pairs – in line with strong E-I connections typically found in the cerebral cortex (Campagnola et al. 2021).

(C) Left: mean peak-aligned EPSCs, pooled from all experiments (14 patched cells) for connections that were classified to be excitatory (n = 75; decay τfast = 2.5 ms, τslow = 15.2 ms) or inhibitory (n = 67; decay τ = 8.6 ms). Right: same as on the left, but individual EPSC waveforms were first normalized with respect to their peak amplitude. Shadings denote the s.e.m.

(D) Fluorescence images of a putative inhibitory and excitatory cell (Z-projections of stitched stack mosaics).

Spiking of individual neurons is dominantly controlled by a few strong incoming connections

(A) Optimization of STP estimates. Bottom: example baseline-subtracted CCGs, based on spike trains of a presynaptic and the postsynaptic unit. Top: EPSC of the corresponding connection. The EPSC onset latency relative to the presynaptic spike time was used to set the start time of the STP quantification window (green), thereby providing a more accurate STP estimate. This approach was particularly important for excitatory connections with variable excEPSC onset latencies; e.g., due to variable axonal length or synapse location within the dendritic tree.

(B) Relationship between STP (optimized estimate) and EPSC amplitude for excitatory (left)/inhibitory (right) connections onto individual postsynaptic cells; lines are linear fits (here curves because of the logarithmic axis for amplitude; data points and linear fit belonging to the same postsynaptic cell are displayed in the same color). Insets show slopes from the linear fits with mean ± s.d. Only cells with at least three E/I inputs were included ([E/I] 5/6 cells and 40/33 connections). Across all connections, the mean STP value for excitatory and inhibitory connections was 0.0025 ± 0.0150 s.d. and -0.0090 ± 0.0099 s.d., respectively.

Data were comprehensively analyzed by linear mixed-effects modeling (with EPSC amplitude as a fixed effect and postsynaptic cell ID as a random effect). Significance was assessed by a likelihood ratio test. *P < 0.05, **P < 0.01

Neuronal spiking is partially governed by rapid and brief changes in synaptic input activity

(A) Top: raster plot of spike times for all presynaptic units of an example postsynaptic cell. Bottom: reconstructed inhibitory (gi; blue) and excitatory (ge; red) synaptic conductance traces. Spike times of the postsynaptic unit in magenta. E/(E + I) ratio (black trace; right).

(B) Mean auto-(ACG) and cross-correlograms (CCG) of the reconstructed giand ge of the neuron in (A) (25 individual mean-subtracted high-conductance events; shadings denote s.e.m.).

(C) Pairwise unit synchrony (STTC with 10 ms binning [Cutts and Eglen 2014]; [E/I] 330/168 unit comparisons; U = 16288, P < 0.001, Mann-Whitney U test). Box plots indicate median and interquartile range, and whiskers the minimum/maximum values except for outliers.

(D) Spike-triggered average of synaptic input conductances and E/(E + I) ratio of an example neuron (n = 15’813 postsynaptic spikes). See also Fig. S3 for variability across spikes and Fig. S4 for contributions of individual inputs.

(E) Mean (black) and individual baseline-subtracted and peak scaled spike-triggered average E/(E + I) traces (n = 7 cells).

(F) E/(E + I)ST A (black; mean from (E)) aligned to intracellular AP (green; mean from Fig. S5). See Fig. S5 for the temporal relationship between extracellularly measured spike times and intracellularly recorded action potential waveforms.

(G/H) Density scatter plots of two example neurons, displaying the changes in input excitation and inhibition that occurred right before individual postsynaptic spikes (for each spike time t, we calculated the mean g value from t-1 ms to t+1 ms and subtracted the mean g value from t-20 ms to t-10 ms). As a reference (see plots labeled ‘non-spiking’), we also assessed the changes in input conductances that occurred during non-spiking periods (using random time points excluding periods ± 20 ms from measured spikes; 10 x the number of spikes).

(I) The percentage of postsynaptic spikes for the four possible pre-spike conductance-change combinations: increase in I and decrease in E (top-left; Z = 2.4, P = 0.018, Wilcoxon signed rank test); increase in I and E (top-right; Z = 0.68, P = 0.50); decrease in I and E (bottom-left; Z = -0.68, P = 0.50); decrease in I and increase in E (bottom-right; Z = -2.4, P = 0.018). Measured spikes (‘spk’) were compared to random time points from non-spiking (‘non-spk’) periods (n = 7 cells).

(J) Violin plots of pre-spike changes in input inhibition (blue) and excitation (red) for individual postsynaptic neurons (n = 7). Distributions of non-spiking periods in grey. In each violin plot, the dark-colored area marks the interquartile range, the circle indicates the median, and the horizontal line marks the mean. All ‘spk’ vs. ‘non-spk’ comparisons significantly different except for excitation of neuron 1, 3, 4; otherwise P < 0.001 (Student’s t-test).

Analyses from (G-J) focused on time periods with synaptic input activity (g > 2 × s.d. of g for both E and I). Across postsynaptic cells, 13’742 ± 7’173 s.d. postsynaptic spikes were recorded, and 15% ± 15% s.d. (min = 4%, max = 45%, median = 7%) occurred during high-activity periods. Some extreme values in (G-H) and (J) are not displayed. *P < 0.05, **P < 0.01, ***P < 0.001.

Network-state-dependent coordination of postsynaptic spiking by inhibitory inputs

(A) For an example postsynaptic cell, spike-triggered average synaptic conductances (ge/i,STA) based on postsynaptic spike times that occurred during periods of either high ([E/I] orange/light-blue) or low (red/blue) input conductance states (mean and s.d.). Right: Magnified region. Same cell as in Fig. 6D.

(B) Mean gi,ST A across experiments (n = 6 postsynaptic cells; shadings denote s.e.m.) for high (light-blue) and low g (blue) states; scaled by their respective negative and positive peak values. Neurons with a pre-spike drop in inhibition were included in this analysis (neuron #1-6 from Fig. 6J).

(C) Comparison of the time difference between the negative gi,ST A peak and postsynaptic spike time (Z = -2.2, P = 0.028, Wilcoxon signed rank test) and the time difference between the negative and positive gi,ST A peaks (Z = 2.2, P = 0.027) for low and high g states.

Recording periods were labeled high g (conductance) state when gi > 3 × s.d. of gi, otherwise low g state. *P < 0.05.

Organization of incoming monosynaptic connections at the level of individual postsynaptic cells

(A) Presynaptic spike rate, amplitude of the EPSC estimate (also converted to conductance) and EPSC onset latency were analyzed.

(B-D) For each property, data were examined at the network level (connections from all 14 paired recordings pooled) and from the perspective of the individual postsynaptic cells. For the network-level analysis (left), a histogram with and without previous log-transformation are shown (black curve: Gaussian fit; arrow head marks the peak). For all properties, including amplitude (U = 1660, P < 0.001, n[E/I] = 75/67), spike rate (U = 601, P < 0.001, n[E/I] = 45/33) and onset latency (U = 739.5, P < 0.001 n[E/I] = 75/67), there was a significant difference between excitatory and inhibitory connections (Mann-Whitney U test). Moreover, all distributions, except for that of inhibitory onset latency (D), were approximately log-normal (Shapiro-Wilk test). For the single-cell-level analysis (right), an example single-cell histogram is shown in addition to a box plot of the skewness s. The skewness of amplitude ([E/I]; Z = 2.8/2.8, P = 0.0051/0.0044, n = 10/11 postsynaptic cells) and spike rate (Z = 2.1/2.7, P = 0.037/0.0076), but not onset latency (Z = 1.4/-1.5, P = 0.17/0.13), was significantly different from zero (Wilcoxon signed rank test).

(E) Relationship between EPSC amplitude and onset delay at single-cell level. Top: scatter plot, with one linear regression fit for each postsynaptic cell. Bottom: slope values of linear fits.

(F) Relationship between EPSC amplitude and presynaptic spike rate at the single-cell level. Amplitude values were log-transformed, otherwise as in (E). Only cells with at least 3 E/I inputs were included (n[E/I] = 10/11 cells).

Data in (E/F) were comprehensively analyzed by linear mixed modeling. Box plots indicate median and interquartile ranges and whiskers the minimum/maximum values except for outliers. *P < 0.05, **P < 0.01

A few key inhibitory hub cells with high spike rates, strong synapses and fast action potential propagation dominate the network

(A) Multiple neurons in the same network were sequentially patched (locations approximated by electrode schematic). Paired patch-clamp and HD-MEA recordings were performed for each cell, followed by identification of their incoming connections using our regression approach. Open triangles and circles mark the positions of excitatory and inhibitory cells, respectively. Overlaid circles/triangles represent neurons that were presynaptic to multiple patched cells. Top inset: virtually identical example extracellular footprints from four separate paired recordings (traces slightly shifted in time for better visibility). Bottom inset: Example EPSCs evoked in different postsynaptic cells by the same excitatory/inhibitory presynaptic cell (marked in the network schematic by single/double asterisk).

(B) Outdegree distribution for the network shown in (A).

(C) Skewness s of the amplitude distribution for the outgoing connections of individual neurons (minimum of 3 connections).

(D) Relationship between the spike rate of a neuron and its outdegree. Significance was assessed after applying the Holm correction for multiple comparisons.

(E) Left: Relationship between the sum of EPSC amplitudes of all outgoing connections of a neuron and its outdegree. Right: for each presynaptic cell, the mean postsynaptic EPSC waveform was computed; the mean of all from cells with outdegree 1-3 and 4-7 are shown.

(F) Left: Relationship between action potential velocity of a presynaptic neuron and its outdegree. Middle: Example HD-MEA action potential (AP) latency plot from a neuron with outdegree 1 and 7. Right: Plots of AP latency vs. distance between electrode and soma for the two example units.

Variations in the number of identified synaptic inputs

Example synaptic input reconstructions of two neurons that were sequentially patched in the same network and yet displayed a difference in the number of identified incoming connections [(A) 20 vs (B) 6]. Example periods of unit spiking activity of the entire network and only of the respective presynaptic neurons (top two raster plots) are shown, in addition to the measured (black) and reconstructed (magenta/green) input-current traces. The current-trace reconstruction showed a good matching with the measured VC patch-clamp recording for both neurons, which suggests that the difference in identified inputs was of biological origin and not due to an incomplete input mapping.

Inhibitory and excitatory neurons show differences in action potential propagation, neuronal foot-print size and axonal reach

(A) Relationship between excEPSC (red) / inhEPSC (blue) onset latency (relative to presynaptic spike time) and the distance between pre- and postsynaptic cell. Lines represent linear regression fits ([E/I] 47/38 connections from 7 cells).

(B) To visualize and quantify action potential (AP) propagation for a given unit, the spike-triggered average for all simultaneously recorded traces (up to 1024) was computed, and the peak latency relative to the electrode featuring the largest signal (marked by a white plus sign) was determined. Two example AP latency plots for inhibitory units with an extensive (top) and relatively small (bottom) axonal arbor are shown. Gray arrow heads in the top plot indicate the electrodes for which the STA traces are displayed underneath (filtered at 0.3–9.5 kHz). Gray dots indicate electrodes that did not show any AP signal. White ’gaps’ indicate electrodes that were not selected for long-term recording.

(C) As in (B), but for excitatory units.

(D) Difference between excitatory and inhibitory units in terms of action potential propagation velocity (U = 340, P < 0.001, Mann-Whitney U test), footprint size (U = 488.5, P = 0.010) and axonal reach (U = 425, P = 0.0013); [E/I] 45/33 pooled presynaptic units; for each unit, 1500 spiking events (# spikes available for all units) were averaged to generate the footprint. AP propagation velocity was quantified based on the AP latencies across the entire unit footprint (see Methods for details). Axonal reach is defined as the 90th percentile of the distances between the electrodes that belong to the unit footprint (i.e., show an AP signature) and the largest signal electrode (presumably near the soma). Box plots indicate median and interquartile ranges, and whiskers the minimum/maximum values except for outliers.

*P < 0.05, **P < 0.01, ***P < 0.001.

Variability in near-spike input conductances across entire recording

The spike-triggered-average conductances, ge/i,STA, from Fig. 6D shown with s.d. (shadings). Note that spike times from both low- and high-activity states were included (n = 15’813 postsynaptic spikes).

Role of individual incoming connections in controlling postsynaptic spiking

The spike-triggered-average conductance, ge/i,STA, was calculated separately for the conductance trace of each incoming connection to an example postsynaptic cell (9 inhibitory and 11 excitatory inputs). There appeared to be a few key inputs that had a particularly strong role in controlling postsynaptic spiking. Moreover, there was a subset of excitatory inputs with a sharp increase in conductance near postsynaptic spiking, while others only contributed through slow conductance changes. Interestingly, postsynaptic spiking typically occurred during the early rising phase of the rapid increases in excitatory conductance. As expected, the contribution of each connection to the total STA conductance depended on the respective synaptic strength and presynaptic spike rate: the mean ge/i,STA (from -2 s to +2 s) for each input correlated well with the respective product of input amplitude and presynaptic spike rate, for both excitation (R2 = 0.77; P < 0.001; n = 11) and inhibition (R2 = 0.96; P < 0.001; n = 9).

Relationship between extracellularly detected unit spike times and corresponding intracellular action potentials

Paired HD-MEA and IC whole-cell patch-clamp recordings during spontaneous spiking activity were performed, and the unit of the patched cell was identified by footprint matching (analogous to Fig. 2D), following spike sorting of the HD-MEA data. An IC internal solution with physiological chloride concentration was used.

(A) Example patch-clamp voltage traces with action potentials. Red triangles mark the peak of the intracellular action potentials. Yellow lines indicate spike times of the unit that corresponded to the patched cell.

(B) Histogram of the time differences between unit spike time and intracellular action potential peak of the same cell as in (A).

(C) Mean action potential waveform and mean unit spike time (shadings: s.d.) of the same cell as in (A/B).

(D/E) Results of two more cells, displayed as in (C).

Overview of changes in input inhibition before and after individual postsynaptic spikes

(A) Analogous to Fig. 6J, the change in input inhibition before individual postsynaptic spikes was quantified for events that occurred during low-activity states. Often narrower ‘spk’ distributions could indicate restricted windows of opportunity for spiking (e.g., only after significant decay of any input inhibition).

(B/C) Change in input inhibition right after postsynaptic spikes that occurred during high- (B) and low-activity (C) states. For this quantification, the mean gi value from t+5 ms to t+10 ms was calculated for each spike time t, followed by subtraction of the mean gi value from t-1 ms to t+1 ms.

(D) Standard deviations across the spike-triggered average input inhibition traces, calculated separately for spiking events from high- (top) and low-activity (bottom) periods. Traces are aligned to the standard deviation value at spike time.

Box plots indicate median and interquartile ranges and whiskers the minimum/maximum values except for outliers (outliers not shown). In the reference distribution of neuron 5 in (A/C), quartiles Q1-Q3 equal zero. By definition, a neuron experienced an inhibitory high-conductance event when gi> 3 × s.d. of gi. For non-spike reference distributions, time points were randomly selected from either high- or low-activity periods, excluding ± 20 ms adjacent to measured spikes. It was n = 6 postsynaptic neurons in all panels.

Spatial distribution of presynaptic unit footprints and monosynaptic connections that were identified based on multiple paired recordings in the same network.

Equivalent to Fig. 9A, with unit footprints shown (filtered at 0.3–9.5 kHz). Only the traces from electrodes near the electrode with the largest signal amplitude are displayed for each unit. Note that some footprints match in terms of both location and signal characteristics, suggesting that the identified units represent the name neuron.

Spike-triggered-average-based estimation of EPSC waveforms using paired HD-MEA and patch-clamp recordings

An alternative to the regression-based EPSC estimation approach relied on calculating, for each unit in the network, the spike-triggered average patch response using only ‘isolated events’; that is, presynaptic spike times around which none of the other potentially connected units exhibited spiking (see Methods for details).

(A) Top: example spike raster plot of presynaptic activity with corresponding whole-cell current-clamp recording. Isolated events are highlighted in colour. Overlapping responses (arrow) were excluded. (‘Unit 0’ contains spike times from units that might be connected to the patched cell, but without reasonable confidence). Bottom: for three example inputs, the mean (colored) and the individual EPSCs (gray) across the entire recording are shown, in addition to the respective EPSC amplitude histogram (a Gaussian filter with a 20-element sliding window was applied to individual response traces for better visualization; mean trace based on raw data).

(B) Another example of evoked EPSCs and amplitude histogram from a separate recording. ‘Peaks’ in the amplitude distribution consistent with multi-site quantal release.

(C) Left: Comparison of EPSC-waveform estimates based on regression vs. STA method of example connections shown in (A). Right: EPSC amplitudes, based on connections identified by both methods, were highly correlated.

Effect of inclusion current thresholds on the estimation of EPSC waveforms

Different current thresholds, that determined which parts of the recording were included in the regression analysis, were applied to the patch-clamp current traces. The threshold values were multiples of the standard deviation of the respective current trace. Specifically, for each recording, 6 threshold values were evaluated: 1, 2, 3, 4, 5 & 10 × 10 × ŝ of the current trace. With ŝ as m.a.d. × 1.4826.

(A) The determined current thresholds across experiments.

(B) We used the EPSCs determined by our STA method, which was performed without any inclusion criteria involving current thresholds, as ground-truth EPSCs. To these ground-truth EPSCs we compared the regression-EPSC estimates that were determined using the different current thresholds. Specifically, for each paired recording, we calculated for each common synaptic input the ratio of the EPSC amplitudes determined by the regression and STA method. The mean of these ratios was calculated for each postsynaptic cell, which is summarized here.

(C) Largest EPSC amplitude for each cell.

(D) Difference between current threshold and largest EPSC amplitude. In all panels, box plots indicate median and interquartile ranges, and whiskers the minimum/maximum values (no outliers). Some extreme whiskers in (A) and (D) were trimmed for better data visibility. [14 cells in (A-D); 9 patched cells in (B), where recordings were only included if at least three common inputs were identified by using the two methods.]

Compared to the EPSC amplitudes determined by the STA method, the different current thresholds had only a small effect on the EPSC amplitudes that were determined by the regression-based method (B). Moreover, the maximum EPSC amplitude that was identified showed little dependence on the current threshold (C). These findings suggested that synaptic inputs exhibited sufficient activity – i.e., for a reliable EPSC estimation – during periods that were selected by even the smallest current thresholds. Moreover, the current threshold values were well above the maximum EPSC amplitudes as confirmed in (D). We chose a threshold of 30 × ŝ of the current trace, which resulted in a median current threshold of approximately 200 pA, with sufficient distance to the maximum EPSC amplitude, while also keeping computational costs within a feasible range.

Estimation of synaptic input waveforms based on current-clamp recordings

Two paired HD-MEA and whole-cell patch-clamp recordings were acquired from the same cell in either voltage-clamp (A) or current-clamp mode (B).

(A/B) Top: raster plot of the identified presynaptic units of an example period (corresponding presynaptic cells are marked by grey arrows). Middle: patch-clamp recording (black) and reconstruction based on synaptic inputs (green). Bottom: location of presynaptic cells; one coloured filled circle per cell (black dots: HD-MEA electrodes selected for recording; grey dots: electrodes not recorded from). The same presynaptic cells identified in (A) and (B) have matching colours. The larger the circle the stronger the connection.

(C) Left: Corresponding EPSP- and EPSC-waveform estimates (waveforms from the same input have matching colors). Right: Correlation of EPSP and EPSC amplitudes (linear regression fit; R2 = 0.99).