Deconvolution improves the detection and quantification of spike transmission gain from spike trains — TEI-tagged view

SpivakpersName, SloinpersName, SomeckpersName, StarkpersName, LevipersName (2022date). “Deconvolution improves the detection and quantification of spike transmission gain from spike trains.” Communications Biologytitle·j.

Object of study — how much of what is read as a connection between two neurons is contributed by the firing pattern of each neuron on its own

Research question — whether EXISTING METHODS FOR DETERMINING NEURONAL CONNECTIVITYterm FROM SPIKE TIMES ARE HIGHLY AFFECTED BY BURSTterm SPIKING ACTIVITY, RESULTING IN OVER- OR UNDERESTIMATION OF SPIKE TRANSMISSIONterm

Contribution — a decomposition of the cross-correlation into three elements and A DECONVOLUTION-BASED ALGORITHM FOR REMOVING EFFECTS OF SECOND-ORDER SPIKE TRAIN STATISTICS

Full text — tagged

The reference list is not part of the text we hold for this paper. Our copy was extracted from the publisher's page and the bibliography did not survive that extraction, so the article below ends before its references. Nothing has been removed by us.

Deconvolutionterm improves the detection and quantification of spike transmissionterm gain from spike trains Lidor SpivakpersName Lidor SpivakpersName 1 Sagol School of Neuroscience and Department of Physiology and Pharmacology, Sackler Faculty of Medicine, Tel Aviv University, Tel Aviv, 6997801 Israel 1 , Amir LevipersName Amir LevipersName 1 Sagol School of Neuroscience and Department of Physiology and Pharmacology, Sackler Faculty of Medicine, Tel Aviv University, Tel Aviv, 6997801 Israel 1 , Hadas E SloinpersName Hadas E SloinpersName 1 Sagol School of Neuroscience and Department of Physiology and Pharmacology, Sackler Faculty of Medicine, Tel Aviv University, Tel Aviv, 6997801 Israel 1 , Shirly SomeckpersName Shirly SomeckpersName 1 Sagol School of Neuroscience and Department of Physiology and Pharmacology, Sackler Faculty of Medicine, Tel Aviv University, Tel Aviv, 6997801 Israel 1 , Eran StarkpersName Eran StarkpersName 1 Sagol School of Neuroscience and Department of Physiology and Pharmacology, Sackler Faculty of Medicine, Tel Aviv University, Tel Aviv, 6997801 Israel 1, ✉ 1 Sagol School of Neuroscience and Department of Physiology and Pharmacology, Sackler Faculty of Medicine, Tel Aviv University, Tel Aviv, 6997801 Israel ✉ Corresponding author.

Received 2022 Jan 20; Accepted 2022 May 4; Collection date 2022.

© The Author(s) 2022 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/ .

PMCID: PMC9156687 PMID: 35641587 Abstract Accurate detection and quantification of spike transmissionterm between neurons is essential for determining neural network mechanisms that govern cognitive functions. Using point process and conductance-based simulationsterm, we found that existing methods for determining neuronal connectivityterm from spike times are highly affected by burstterm spiking activity, resulting in over- or underestimation of spike transmissionterm. To improve performance, we developed a mathematical framework for decomposing the cross-correlation between two spike trains. We then devised a deconvolution-based algorithm for removing effects of second-order spike train statistics. Deconvolutionterm removed the effect of burstterm spiking, improving the estimation of neuronal connectivityterm yielded by state-of-the-art methods. Application of deconvolutionterm to neuronal data recorded from hippocampal region CA1term of freely-moving mice produced higher estimates of spike transmissionterm, in particular when spike trains exhibited burststerm. Deconvolutionterm facilitates the precise construction of complex connectivityterm maps, opening the door to enhanced understanding of the neural mechanisms underlying brain function.

Subject terms: Neural circuits, Neural decoding, Neural decoding A mathematical framework and algorithm are developed that decompose the cross-correlation histogram between two spike trains to determine connectivityterm between two neurons.

Introduction Neurons in the brain act as information gates, each receiving inputs from thousands of different sources and deciding which to pass on to other neurons. The capability of transmitting relevant inputs to other cells is at the very core of sensory processing, action, and cognition. One way for investigating how information is transmitted between two neurons is to measure the postsynaptic potential (PSP) in one cell following a spike of a presynaptic cell 1 – 3 . When excitatory PSP (EPSPterm) magnitude is larger, the probability that the postsynaptic neuron will fire following a presynaptic spike is higher 2 . However, EPSPterm magnitude is affected not only by synaptic strength, but also by the instantaneous membrane potential of the postsynaptic neuron 4 . Moreover, even constant EPSPterm magnitude cannot predict whether a postsynaptic neuron will fire after a presynaptic spike. When an EPSPterm is generated at lower membrane potentials, cells are less likely to fire 5 . Thus, EPSPterm magnitude alone is insufficient for determining whether a postsynaptic neuron will fire following a presynaptic spike, thereby allowing information to pass on to other cells.

An alternative approach for determining connectivityterm between two neurons is to directly quantify whether presynaptic spikes generate spikes in the postsynaptic neuron. The spike-to-spike cross-correlation histogram (CCHterm) 6 – 8 has been used for measuring the transmission of spikes between two neurons 9 – 11 . However, most CCHterm applications employ spontaneous spiking, and are therefore correlative measurements which cannot differentiate causal transmission from non-causal spike patterns. For example, firing of two cells may be co-modulated by another source such as visual input or a third neuron, causing the two cells to fire in synchrony 12 . Thus, the CCHterm is influenced not only by spike transmissionterm properties but also by co-modulation of the pre- and postsynaptic trains 13 – 15 .

To differentiate between direct spike transmissionterm due to inter-neuronal coupling and other indirect sources, several methods have been proposed 16 – 20 . An underlying assumption in these approaches is that while direct connections exhibit fast transients in the CCHterm, indirect sources produce slower fluctuations. Based on this assumption, timescale separationterm can be used to detect and estimate synaptic connectivityterm, differentiating between fast transients (which are presumably due to direct connectivityterm), and slower fluctuations (which are presumably due to other sources). Indeed, when timescale separationterm methods are employed, spike transmissionterm measurements derived from casual methods are consistent with measures derived from spontaneous spiking CCHsterm 21 . However, it is still unclear to what extent timescale separationterm methods can differentiate between transmitted spikes and third-party, synchronous spiking. Moreover, the CCHterm is influenced not only by the interaction between the pre- and postsynaptic spike trains, but also by the spiking activity pattern of each neuron 8 . When one of two connected neurons exhibits high-frequency periodicityterm or burstterm spiking, fast features appear in the CCHterm 22 , 23 . These features cannot be differentiated from transmitted spikes by timescale separationterm. Thus, even when timescale separationterm methods are used, transmitted spikes cannot be fully differentiated from non-transmission patterns such as burstterm spiking activity.

To remove non-transmission features from the CCHterm and improve the estimation of functional connectivityterm from spike trains, we formulated a mathematical framework for CCHterm composition. We found that the CCHterm can be expressed as a sum of three distinct elements. Based on this finding, we developed a deconvolution-based algorithm for removing second-order spike train statistics from the CCHterm. Using deconvolutionterm in concert with timescale separationterm methods removed the effect of burstterm spiking, improving the accuracy of both detection and quantification of neuronal connectivityterm in synthetic data. Consistent with simulationterm results, application to neuronal data recorded from CA1term of freely-moving mice showed that deconvolutionterm increased estimates of spike transmissionterm, in particular when the spike trains exhibited burststerm.

Results Spike transmissionterm gain between two neurons can be accurately estimated from spike-to-spike cross-correlation Neurons can transmit spiking information to other neurons via different connections including chemical synapses (e.g., AMPA, GABA, NMDA, or other receptors) 24 , 25 , electrical synapses (gap junctions) 26 , or ephaptic coupling (electromagnetic fields) 27 . Whether spiking information will propagate in the brain depends on the effective connectivityterm between neurons. When membrane potentials are accessible, connectivityterm between two neurons can be quantified by the magnitude of the EPSPterm in the postsynaptic neuron following a single presynaptic spike 1 – 3 . However, EPSPterm magnitude does not indicate whether the postsynaptic neuron will pass information onward to other neurons in the form of spikes. Alternatively, the spike trains themselves may be used for estimating the effective connectivityterm between two neurons, defined here as spike transmissionterm.

To determine the relations between EPSP-based estimates and spike transmissionterm, we simulated a two-cell network in which coupling strength was varied (Fig. 1 ). The presynaptic neuron exhibited Poisson spiking modified by a 2 ms refractory period ( λ 1 = 5 spk/s). The postsynaptic neuron was modeled as a conductance-based leaky integrate and fire (LIF) neuron. In the lack of connectivityterm, the postsynaptic neuron did not emit any spikes when membrane potential variability, quantified using σ , was low (Fig. 1a ). When σ was increased, spontaneous spike rate increased monotonically (Fig. 1a ). Specifically, in the range of noiseterm observed in intact preparations ( σ < 4 mV) 5 , 28 , spike rate was in the range observed in freely-moving animals (<25 spk/s) 29 .

Fig. 1. Spike transmissionterm between two neurons can be estimated from spike times.

a Increased membrane potential variability increases spontaneous spiking. A leaky integrate and fire (LIF) model neuron was simulated with different levels of membrane potential variability, quantified by the noiseterm SD, σ ( n = 6 repetitions for each of n = 11 σ values, from 0 to 5 mV with 0.5 mV increments; each 180 min). The firing rate ( λ 2 ) of the simulated neuron is shown for each σ value; here and in ( b ), error bands indicate SEM. The trace next to the cartoon shows 300 ms of simulated membrane potential, and the histograms show the distribution of membrane potentials ( σ = 3 mV). b Even small-magnitude EPSPsterm induce observable spike transmissionterm gain in the presence of noiseterm. The simulated LIF neuron (2) was driven via an excitatory synapse by a presynaptic neuron (1) with Poisson spiking ( λ 1 = 5 spk/s) modified by refractoriness. Synaptic conductance was modified to yield n = 11 different unitary EPSPterm (uEPSP) magnitudes (from 0 to 2 mV with 0.2 mV increments), and the real spike transmission gain (rSTG)term·def was computed for every σ and uEPSP combination ( n = 121). Both uEPSP magnitude and membrane potential variability control spike transmissionterm. c A cross-correlation histogram (CCHterm; bins size: 1 ms) was computed from a pair of simulated spike trains as in ( b ); σ = 3.4 mV, uEPSP magnitude = 0.6 mV. For deriving an estimated STGterm (eSTGterm) from the CCHterm, the level of baseline joint counts (blue line) was first estimated using the “tailsterm” predictor, averaging CCHterm values at time lags | τ | ≥ 11 ms. The spike transmissionterm curve (STC; black line) was estimated as the above-baseline curve with a peak in the temporal region of interest (ROI; defined here as 0 < τ ≤ 5 ms; gray bins), with zeros in all other bins. The eSTGterm is the area under the estimated STC. d Spike transmissionterm gain estimated from the CCHterm is an accurate estimate of the real STGterm. For every simulationterm run ( n = 6 repetitions for each of 121 σ and uEPSP combinations), rSTGsterm and eSTGsterm were computed. Different colors correspond to different uEPSP magnitudes as in ( b ), and the black circle corresponds to the example in ( c ). Spearman’s correlation coefficient between eSTGterm and rSTGterm is ρ = 0.94 ( n = 726 random samples, p < 0.001, permutation testterm·stat).

Next, we connected the pre- and postsynaptic neurons via a synaptic conductance. Conductance magnitude was calibrated to produce unitary EPSPsterm (uEPSPs) in the range of 0–2 mV, as observed in vivo 2 , 28 , 30 , 31 . When membrane potential variability was low ( σ < 2 mV), even the stronger uEPSPs failed to generate postsynaptic spiking. However, when σ was higher (2–4 mV), uEPSPs of the same magnitude induced postsynaptic spikes in addition to the spontaneous noise-induced spikes (Fig. 1b , red curve). To quantify effective connectivityterm, we defined “real spike transmissionterm gain” (rSTGterm) as the number of extra postsynaptic spikes generated following every presynaptic spike. In the range of biologically-relevant noiseterm, rSTGterm increased monotonically with σ . Furthermore, for every σ , rSTGterm increased monotonically with uEPSP magnitude (Fig. 1b , top inset). While even relatively-large EPSPsterm fail to produce postsynaptic spikes in the lack of membrane potential variability, weaker EPSPsterm (0.2–1 mV) could induce consistent spiking when noiseterm was sufficiently high. Moreover, the uEPSP magnitude required to generate a fixed rSTGterm decreased monotonically with σ (Fig. 1b , bottom inset). Thus, EPSPterm magnitude alone does not suffice for estimating spike transmissionterm.

To determine whether effective connectivityterm associated with weak EPSPsterm can be detected from spike timing alone, we constructed CCHsterm for pre- and postsynaptic spike train pairs (Fig. 1c ). We defined the “spike transmissionterm curve” (STC) as the impulse response of spike transmissionterm between the pre- and postsynaptic neurons. Empirically, an estimated STC (eSTC) can be determined by the difference between the CCHterm and the baseline. Here, we defined the baseline using the “tailsterm” predictor, namely the CCHterm count at longer time lags (| τ | ≥ 11 ms). We then defined the estimated STGterm (eSTGterm) as the area under the STC peak in a temporal region of interest (ROI; defined here as 0 < τ ≤ 5 ms 32 ; Fig. 1c , gray bins), divided by the number of presynaptic spikes. For an STC below baseline, the STGterm is a negative quantity. Previous studies have referred to the STGterm using various other names: “effectiveness” 33 , “asynchronous gain” 34 , “synaptic efficacy” 11 , and “spike transmissionterm probability” 21 , 32 . For the CCHterm in Fig. 1c , generated using relatively weak EPSPsterm (0.6 mV; σ = 3.4 mV), the eSTGterm (0.032) was of the same order of magnitude as the rSTGterm (0.028, 116%; Fig. 1c ). Considering all non-zero rSTGterm values, the mean eSTGterm error (eSTGterm–rSTGterm)/rSTGterm, was 2.5% (456 spike train pairs with non-zero rSTGterm; smaller than 5%, p < 0.05, one-tailed Wilcoxon’s signed-rank test; Fig. 1d ). Thus, the STGterm estimated from the CCHterm between spike trains closely matches the real spike transmissionterm gain.

The cross-correlation between two spike trains can be expressed as a sum of three elements The foregoing discussion (and Fig. 1 ) focused on a specific configuration in which the presynaptic spike train was generated via a homogenous Poisson process modified by refractoriness. Furthermore, all statistical dependencies between the two spike trains were captured by a single modeled synapse. In general, dependencies within and between spike trains may be richer (Fig. 2a ). Individual spike trains typically deviate from Poisson processes, conforming to higher-order gamma processes or exhibiting spike burststerm 35 , 36 . Connectivityterm between neurons is often not purely feedforward, exhibiting reciprocal excitation 1 or feedback inhibition 37 , 38 . Furthermore, dependencies between two spike trains may result from covariance of other inputs to the two neurons 6 , 8 , 14 , 39 rather than from synaptic connections between the neurons. These and other possible sources of complexity give rise to the observed spike trains (Fig. 2b ) from which the CCHterm is derived (Fig. 2c ). Therefore, the STC estimated from a pair of spike trains is generally not identical to the real STC. As a direct consequence, the eSTGterm is not identical to the rSTGterm.

Fig. 2. The cross-correlation between two spike trains can be expressed as a sum of three elements.

a Cartoon of two neurons, driven by external sources and coupled by reciprocal excitatory and inhibitory monosynaptic connections. In the example, both neurons exhibit spiking activity which deviates from a Poisson process: the excitatory neuron exhibits burstterm spiking, and the inhibitory neuron exhibits second-order gamma spiking activity. b The spike trains of the two neurons result from all of these (and possibly other) sources. c An example CCHterm created from the synthetic spike trains described in ( a , b ). The rSTGterm is 0.04, whereas the eSTGterm provides an inaccurate estimate (0.054; 135%). The example demonstrates that in the presence of complex dynamics, the STGterm cannot be estimated accurately by the “tailsterm” method. d The CCHterm between two coupled spike trains, s 1 and s 2 , can be expressed as a sum of three elements. The first element is the convolutionterm of the auto-correlation histogram of the first train (ACH 1 ) with the STC from neuron 1 to neuron 2 (STC 21 ). The second element, denoted as “other inputs”, is equal to the cross-correlation between the uncoupled spike trains. The third element is the convolutionterm of ACH 2 with STC 12 . e Based on the composition ( d ), a decomposition procedure may recover the other inputs and the two STCs from the CCHterm.

In general, the STC is not necessarily identical to the difference between the CCHterm and the baseline. However, if the STC is fixed (more precisely, if the system is linear time-invariant), there is an exact mathematical dependency between the STC and the CCHterm. In the absence of external input to the postsynaptic neuron, all postsynaptic spikes ( s 2 ) are due to spontaneous presynaptic spikes ( s 1 0 ). In that case, the CCHterm is identical to the convolutionterm (*) between the spike-to-spike auto-correlation histogram (ACH) of neuron 1 spikes, and the STC from neuron 1 to neuron 2 (Fig. 2d ): CCHterm = ACH 1 * STC 21 1 When external input is added to the postsynaptic neuron ( s 2 0 ), the CCHterm is equal to the sum of the above convolutionterm and the cross-correlation (⋆) between the spike trains which are not due to the direct connectivityterm: CCHterm = ACH 1 * STC 21 + s 1 0 ⋆ s 2 0 2 When the connectivityterm between the neurons is bidirectional, and external (possibly correlated) inputs are added to both neurons, the CCHterm can be approximated as a sum of three elements: CCHterm = ACH 1 * STC 21 + s 1 0 ⋆ s 2 0 + ACH 2 * STC 12 3 The approximation is valid as long as the product of the two STGsterm (i.e., the loop gain) is small. Thus, a CCHterm can be described as a sum of three elements: (1) the convolutionterm between the ACH of neuron 1 and the impulse response from neuron 1 to neuron 2; (2) the convolutionterm between the ACH of neuron 2 and the impulse response from neuron 2 to neuron 1; and (3) the cross-correlation between the background spike trains, those not due to direct connectivityterm between the two neurons. The last expression has the conceptual benefit of opening the loop, effectively describing a feedback system (Fig. 2a ) as a feedforward system using an arithmetic decomposition (Fig. 2e ). Moreover, understanding that the CCHterm is a sum of three elements allows reconstructing hidden elements, namely the STCs, from the CCHterm using numerical methods.

Deconvolutionterm eliminates burstterm spiking effects and enables accurate STGterm and PSP estimation To estimate inter-neuronal effective connectivityterm from spike trains, several methods have been proposed 20 , 32 , 40 – 42 . These and other procedures rely on the concept of timescale separationterm, in which processes that occur at fast timescales (e.g., monosynaptic spike transmissionterm) are distinguished from processes occurring at slower timescales (e.g., co-modulated other inputs, such as external inputs or network oscillations). However, in addition to other inputs, the CCHterm is comprised of the convolutionterm products of each ACH with the corresponding STC (Fig. 2d ), which occur at fast timescales. Hence, even if an optimal estimate of the other inputs is removed, a CCHterm between two synaptically-connected neurons will still contain contributions from the ACHs of each neuron.

To demonstrate the effect of second-order spike train statistics on connectivityterm estimation, we simulated two spike trains coupled by an excitatory monosynaptic connection (rSTGterm = 0.04; 833 min; Fig. 3a–c ). The spike train of the presynaptic neuron was initially simulated as a Poisson process modified by an absolute refractory period (ARP; ARP 1 = 2 ms; λ 1 = 2 spk/s), and exhibited bursting activity (Fig. 3b ). The extent of bursting was controlled by a “burst fractionterm·def” (BR) parameter, defined as the probability of each spike to be followed by another spike within a short inter-spike interval (3–7 ms). BR was set to 0.4. The spike train of the postsynaptic neuron realized a second-order gamma process modified by refractoriness ( λ 2 = 8 spk/s; ARP 2 = 2 ms; Fig. 3a ). Consistent with the construction, the CCHterm between the coupled spike trains exhibited a peak within the causal ROI (Fig. 3c ). In addition, the CCHterm exhibited side lobes: secondary peaks with increased counts on both sides of the causal peak, which are due to the bursting activity of the presynaptic neuron. Similar patterns were previously reported in real datasets 22 .

Fig. 3. In the presence of spike burststerm, deconvolutionterm enables accurate recovery of spike transmissionterm gain.

a – e Simulated spike trains of two neurons coupled by an excitatory monosynaptic connection. a ACH of the postsynaptic neuron. b ACH of the presynaptic neuron, exhibiting bursting activity. “Burst fractionterm·def” (BR; the probability of each spike to be first in a burstterm) is 0.4. c CCHterm between the two neurons. In addition to a peak within the causal ROI (0< τ ≤5 ms), the CCHterm exhibits side lobes which are due to the burstterm spiking of the presynaptic neuron. d Conditional rate histograms, derived from the CCHterm in ( b ). e CCHterm with slow (black) and fast GLM coupling filters (pink and gray). f Burstterm spiking affects the ACH ( b ) and the CCHterm ( c ) and causes inaccurate STGterm estimates ( d ). The scheme illustrates the deconvolutionterm process for removing ACH artifacts from the raw CCHterm. g A deconvolved CCHterm·def (dcCCH; green) is derived from the ACHs ( a , b ) and CCHterm ( c ). The dcCCH is free from the effect of spike burststerm, recovering the unidirectional spike transmissionterm curve STC 21 used in the simulationterm. h Conditional rate histograms. i dcCCH with PSPs estimated by the CoNNECT and GLMCCterm methods. j , k CCHterm deconvolutionterm improves STGterm estimation in the presence of burstterm spiking. Spike trains coupled by an excitatory monosynaptic connection were simulated while modifying BR ( n = 30 repetitions at n = 9 BR values, from 0 to 0.4 with 0.05 increments). j eSTG-to-rSTG ratio as a function of BR for the three STGterm estimation methods using CCHsterm (solid lines) and using dcCCHs (dashed lines). In the presence of bursting, deconvolutionterm improves STGterm estimation. k Estimated PSP as a function of burstiness for the CoNNECT and GLMCCterm methods using CCHsterm and dcCCHs. l Normalized slope of the best linear fitterm·stat of the curves in ( j , k ). Error bars, SEM. With the exception of CoNNECT, deconvolutionterm reduces the effect of bursting for all methods ( n = 270 random samples; n.s./*** p > 0.05/ p < 0.001, bootstrap test, n = 300 iterations). m – o Spike trains were coupled by an inhibitory monosynaptic connection. m Example CCHterm. n dcCCH for the same example. o dcCCH with slow GLM (black) and fast GLM coupling filters (pink and gray). p , q Same as ( j , k ), for simulated inhibitory connections ( n = 30 repetitions at every BR value). r Same as ( l ), for inhibitory connections ( n = 270 random samples; n = 300 bootstrap iterations). s – v Two spike trains with burststerm were coupled by an excitatory monosynaptic connection. s ACH of the postsynaptic neuron (BR = 0.1). t ACH of postsynaptic neuron (BR = 0.4). u CCHterm of the two neurons. v dcCCH for the CCHterm in ( u ). w , x Same as ( j , k ), for two spike trains with burststerm ( n = 30 repetitions at every BR value). y Same as ( l ), for two spike trains with burststerm ( n = 270 random samples; n = 300 bootstrap iterations).

To quantify the effect of the ACH on the estimation of effective connectivityterm, we used five different methods: tailsterm, jitterterm, medianterm, GLMCCterm 40 , and CoNNECT 41 . While the tailsterm, jitterterm, and medianterm filtering methods yield estimates of effective connectivityterm measured as gain (STGterm; Fig. 1 ), the GLMCCterm and CoNNECT methods produce estimates in units of mV (PSP). The estimates yielded by the PSP-based methods cannot be gauged directly in the case of the synthetic CCHterm depicted in Fig. 3c , since no ground truthterm for PSP exists in the point process simulationsterm. The tailsterm predictor overestimated the rSTGterm (eSTGterm = 0.053, 133%; Fig. 3d ), and the jitterterm and medianterm predictors underestimated the rSTGterm (jitterterm: 0.023, 58%; medianterm filtering: 0.028, 70%). Thus, when burstterm spiking activity is present, the tailsterm, jitterterm, and medianterm methods fail to produce an accurate eSTGterm.

For accurate estimation of effective connectivityterm from the CCHterm in the presence of arbitrary second-order spike train statistics (e.g., bursting or periodicityterm), we developed a deconvolutionterm algorithm that takes advantage of the CCHterm mathematical properties described in Fig. 2d, e . In a nutshell, the algorithm mitigates the effect of bursting activity by deconvolving the ACHs from the CCHterm, resulting in a CCHterm which does not exhibit ACH-induced side lobes (Fig. 3f ). The outcome is a “deconvolved CCHterm·def” (dcCCH). The dcCCH still includes the additive effect of other inputs (Fig. 2d ). However, the other inputs element can be estimated using timescale separationterm methods and removed by subtraction. Thus, deconvolving the ACHs from the CCHterm before applying timescale separationterm removes ACH artifacts, and may improve connectivityterm estimates.

To test whether deconvolutionterm indeed improves STGterm/PSP estimation, we first applied the algorithm to the CCHterm and ACHs of Fig. 3a–c and re-estimated the STGterm and PSP from the resulting dcCCH (Fig. 3g ). Compared to the estimates derived from the CCHterm (Fig. 3d, e ), all STGterm estimates based on the dcCCH were improved (Fig. 3h ), and both PSP estimates were higher (Fig. 3i ). Specifically, the tailsterm predictor yielded an eSTGterm of 0.0408 (102% of the rSTGterm); the jitterterm predictor underestimated the rSTGterm (0.0276, 69%); and medianterm filtering yielded an eSTGterm of 0.0392 (98%; Fig. 3h ). We then carried out simulationsterm in which the presynaptic burst fractionterm·var was varied in an orderly manner (from 0 to 0.4 with 0.05 increments; rSTGterm = 0.04; 30 repetitions each; other parameters were kept the same as in Fig. 3a–c ). For each simulationterm, estimated STGsterm and PSPs were computed from the CCHterm and dcCCH using all five methods (Fig. 3j, k ). For the tailsterm, jitterterm, and medianterm predictors, eSTGsterm derived from the dcCCH were more accurate than the eSTGsterm derived from the CCHterm when burst fractionsterm·def were above 0.15 ( p < 0.001, Mann–Whitney U -test). For the GLMCCterm method, only the PSP estimations derived from the raw CCHterm depended on the burst fractionterm·def (GLMCCterm without deconvolutionterm: p < 0.001; with deconvolutionterm: p = 0.13; permutation testterm·stat). For comparing STGterm and PSP-based estimates, we used the normalized slope of the best linear fitterm·stat, after dividing all values by the mean estimates at zero burst fractionterm·def. With the exception of CoNNECT, deconvolutionterm reduced the effect of bursting for all methods ( p < 0.001, bootstrap test; Fig. 3l ). Thus, in the presence of burstterm spiking activity, deconvolutionterm improves STGterm and PSP estimation of excitatory monosynaptic connections.

Second, to test whether deconvolutionterm also improves the estimation of inhibitory monosynaptic connections, we repeated the simulationsterm of Fig. 3g–l with inhibitory connections (rSTGterm = −0.02; Fig. 3m–r ). The spike train of the presynaptic neuron realized a second-order gamma process modified by refractoriness ( λ 1 = 8 spk/s; ARP 1 = 2 ms). The spike train of the postsynaptic neuron was simulated as a Poisson process modified by refractoriness ( λ 2 = 2 spk/s; ARP 2 = 2 ms) which also exhibited bursting activity, varied in an orderly manner. For the tailsterm, jitterterm, and medianterm predictors, eSTGsterm derived from the dcCCH were more accurate than the eSTGsterm derived from the CCHterm when burst fractionsterm·def were above 0.15 (Fig. 3p ; p < 0.001, U -test). We then computed the normalized slope of the best linear fitsterm·stat (Fig. 3r ). For all methods except CoNNECT, deconvolutionterm reduced the effect of bursting (tailsterm, jitterterm, medianterm: p < 0.001; GLMCCterm: p = 0.015; bootstrap test). Thus, in the presence of burstterm spiking activity, deconvolutionterm improves STGterm and PSP estimation for inhibitory monosynaptic connections.

Third, to test whether deconvolutionterm improves the estimation in the case of two directly-connected neurons which both exhibit bursty behavior, we simulated two bursting neurons coupled by a monosynaptic excitatory connection (rSTGterm = 0.04; Fig. 3s–u ). The burst fractionterm·def of the postsynaptic neuron was 0.1 (Fig. 3s ), and the burst fractionterm·def of the presynaptic neuron was varied from 0 to 0.4 with 0.05 increments. For every simulationterm, we estimated STGsterm and PSPs from the CCHterm and dcCCH using all five methods (Fig. 3u–x ). We found that for the tailsterm, jitterterm, and medianterm predictors, eSTGsterm derived from the dcCCH were more accurate than eSTGsterm derived from the CCHterm for burst fractionsterm·def above 0.1 ( p < 0.001, U -test). For the CoNNECT method, only the ePSPsterm derived from the raw CCHterm depended on the burst fractionterm·def (without deconvolutionterm: p = 0.01; with deconvolutionterm: p = 0.053; permutation testterm·stat). For comparing STGterm and PSP-based estimates, we used the normalized slope of the best linear fitterm·stat, after dividing all values by the mean estimates at zero burst fractionterm·def. Deconvolutionterm reduced the effect of bursting for all methods (CoNNECT: p = 0.025; all other methods: p < 0.001, bootstrap test; Fig. 3y ). Thus, when both neurons exhibit burstterm spiking activity, the deconvolutionterm improves STGterm and PSP estimation.

Deconvolutionterm improves effective connectivityterm estimation and detection To evaluate the contribution of deconvolutionterm to the performance of the various methods in a more realistic, noisy scenario, we generated 1250 individual neuronal pairs. Connectivityterm strengths, burst fractionsterm·def, co-modulation, and recording durations were varied widely between pairs (Fig. 4a ). A thousand pairs were generated with either an excitatory or inhibitory connection, and 250 pairs were unconnected. For the excitatory and inhibitory pairs, connection strengths were drawn randomly from a log-normal distributionterm·stat 43 – 45 (mean ± SD excitatory rSTGterm: 0.019 ± 0.01; inhibitory: −0.014 ± 0.007; Fig. 4b ). Burst fractionsterm·def were drawn randomly from a uniform distribution in the [0 0.4] range. Recording duration was set randomly between 90 min and five hours.

Fig. 4. Deconvolutionterm improves quantification and detection of inter-neuronal connections.

a , b A thousand pairs of neurons with monosynaptic connections and 250 unconnected pairs were simulated. Connectivityterm strengths (rSTGterm), burst fractionsterm·def, co-modulation, and recording durations were varied between pairs. a Example CCHsterm from the noisy point process simulationsterm, sorted by rSTGterm values. BR, burst fractionterm·def. b Empirical distribution of the n = 500 excitatory (red) and n = 500 inhibitory (blue) rSTGsterm of the pairs used in the simulationsterm. c – g STGsterm or PSPs are plotted against the ground truthterm (rSTGterm) for dcCCHs ( n = 1000 connected pairs), estimated using the tailsterm ( c ), jitterterm ( d ), medianterm ( e ), GLMCCterm ( f ), and CoNNECT ( g ) methods. Colored (or gray) dots depict estimated STGsterm or PSPs which were (or were not) detected as connections. In all methods (except CoNNECT), the significance threshold used was α = 0.001. In ( c – e ), solid black lines represent perfect estimation (i.e., eSTGterm = rSTGterm). h Mean square errorsterm·stat (MSEterm) for the STG-based methods, in which the ground truthterm and the estimates have the same units. Error bars, SEM. Deconvolutionterm reduces the MSEterm for all STG-based methods ( n = 1000 connected pairs; *** p < 0.001, Wilcoxon’s paired signed-rank testterm·stat). The smallest errors are obtained for medianterm filtering with deconvolutionterm. i – m STGterm ( i – k ; or PSP, l , m ) ratio as a function of burst fractionterm·def. STGterm ratio is defined as the STGterm estimated from the CCHterm, divided by the STGterm estimated from the dcCCH. The PSP ratio is defined in an analogous manner. For all methods except CoNNECT, deconvolutionterm reduces the effect of burstiness. n Slope of best linear fit, computed for ( i – m ) ( n = 1000 connected pairs; n.s./*** p > 0.05/ p < 0.001, bootstrap test, n = 300 iterations). o – t Effect of deconvolutionterm on detection. o – s False-positive (FP) and false-negative (FN) rates were computed for excitatory (Exc.) and inhibitory (Inh.) connections, using CCHsterm (dark bars) and dcCCHs (light bars). Fractions are out of a population of {Exc. FP: n = 2000; Exc. FN: n = 500; Inh. FP: n = 2000; Inh. FN: n = 500} connections. t Mean f 1 scores for all connections, with and without deconvolutionterm. Error bars, SEM. Deconvolutionterm improves detection performance for the tailsterm, medianterm, and GLMCCterm methods (n.s./*** p > 0.05/ p < 0.001, bootstrap test, n = 3500 iterations).

First, we compared estimation performance with and without deconvolutionterm by computing the mean square error (MSE)term·stat for the STG-based methods, for which a ground truthterm exists (Fig. 4c–e ). We found that deconvolutionterm reduced the MSEterm for all STG-based methods ( p < 0.001; Wilcoxon’s paired signed-rank testterm·stat; Fig. 4h ). Of the three STG-based methods, medianterm filtering with deconvolutionterm had the lowest MSEterm (1.27 × 10 −5 ). Medianterm without deconvolutionterm and jitterterm with deconvolutionterm had the second and third lowest MSEsterm, respectively (medianterm: 2.73 × 10 −5 ; jitterterm with deconvolutionterm: 3.48 × 10 −5 ; jitterterm: 5.34 × 10 −5 ). Thus, deconvolutionterm reduced the errors approximately two-fold. To evaluate the effect of deconvolutionterm for all five methods, we computed the estimated STGterm (or PSP) ratio for both excitatory and inhibitory connections. For the “tailsterm” method, the eSTGterm ratio increased with burst fractionterm·def ( ρ = 0.297, p < 0.001, permutation testterm·stat; Fig. 4i ). For the jitterterm, medianterm, and GLMCCterm methods, the estimated STGterm (or PSP) ratio decreased as burst fractionterm·def increased ( p < 0.001, permutation testterm·stat; Fig. 4j–l ). For all methods (except CoNNECT), the slope of the estimated STGterm (or PSP) ratio was different from zero ( p < 0.001, bootstrap test; Fig. 4n ). Together, the reduction of the MSEterm (Fig. 4h ) and the correlation between estimated STGterm (or PSP) ratio and burst fractionterm·def (Fig. 4n ) show that deconvolutionterm reduces the effect of burstterm spiking and improves connectivityterm estimation.

Second, we compared detection performance for all methods (Fig. 4o–t ). We used an f 1 score that combines the effect of false-positive (FP) and false-negative (FN) rates across multiple populations (excitatory and inhibitory connections; Fig. 4t ). Both the medianterm and GLMCCterm, when employed with α = 0.001, yielded low FP rates for excitatory and inhibitory connections. For medianterm filtering, the FP rate for excitatory connections with/without deconvolutionterm was 0.002/0.0012; for inhibitory connections, the FP rates (with and without deconvolutionterm) were 0.0004 (Fig. 4q ). For the GLMCCterm method, the FP rates for the excitatory connections with and without deconvolutionterm were zero; for inhibitory connections, deconvolutionterm improved the FP rate from 0.0004 to zero (Fig. 4r ). Overall, deconvolutionterm improved detection performance for the tailsterm, medianterm, and GLMCCterm methods ( p < 0.001, bootstrap test; Fig. 4t ). Of all methods, GLMCCterm with deconvolutionterm yielded the highest f 1 score (mean [SEM]: 0.955 [0.004]). GLMCCterm without deconvolutionterm and medianterm filtering with deconvolutionterm yielded the second and third scores (GLMCCterm: 0.94 [0.005]; medianterm with deconvolutionterm: 0.92 [0.006]; medianterm: 0.90 [0.007]). Thus, deconvolutionterm improves connectivityterm estimation and detection even for noisy CCHsterm.

Deconvolutionterm improves connection quantification and detection in networks of conductance-based CA1term model neurons The noisy point process simulationsterm (Fig. 4 ) contained well-defined ground truthterm rSTGsterm, and generated CCHsterm that exhibited complex features and multiple types of spike-to-spike interactions. However, the point process simulationsterm were limited to two-neuron “networks” with feedforward connectivityterm and may not exhibit other features, produced in larger-scale networks with reciprocally-connected neurons. To examine the contribution of deconvolutionterm in such settings, we conducted simulationsterm of noisy networks of conductance-based CA1term model neurons. In these simulationsterm, the excitatory neurons (E-cells) were modeled using spiking models that allow incorporating graded bursting behavior (Fig. 5a ). Bursting was controlled by a “burstterm factor” (BF) parameter, set to be between −42 mV (lower bursting activity) and −35 mV (higher bursting activity). Inhibitory neurons (I-cells) did not exhibit spontaneous bursting. In each simulationterm, a hundred-neuron network was constructed by connecting 80 E-cells and 20 I-cells. Connectivityterm patterns were implemented as often assumed for CA1term, with E-to-I, I-to-E, and I-to-I connections, but without E-to-E connections 46 . The E- and I-cell ACHs produced under these conditions (Fig. 5b, c ) are similar to those observed for pyramidal cells (PYR) and interneurons (INT) in real CA1term data (Fig. 6 ).

Fig. 5. Deconvolutionterm improves connection quantification and detection in networks of conductance-based CA1term model neurons.

a Traces show example membrane potentials of six simultaneously-simulated cells. Gray, inhibitory neurons (I-cells); orange, excitatory neurons (E-cells), sorted according to burstiness, quantified by the burstterm factor (BF). b , c Example ACHs and CCHsterm for two E-to-I cell pairs. b A pair with feedforward excitation (Gie = 0.055 mS/cm 2 , Gei = 0 mS/cm 2 ). c A pair with reciprocal excitation and inhibition (Gie = 0.055 mS/cm 2 , Gei = 0.25 mS/cm 2 ). Both pairs were embedded in the same hundred-cell network; in both cases, deconvolutionterm increased the eSTGsterm. d – h A network of 100 neurons (80/20 E/I cells) with fixed connectivityterm strengths was simulated. E-cell BFs were drawn uniformly from the [−42 −35] mV range. d Hinton diagram of network connectivityterm. Red/blue squares represent excitatory/inhibitory connections, and square size is proportional to connectivityterm strength. Dashed lines separate E- and I-cells. e , f eSTGsterm of n = 80 fixed E-to-I connections vs. BF using medianterm filtering without ( e ) and with ( f ) deconvolutionterm. Only eSTGsterm estimated without deconvolutionterm are correlated with the BF. g , h eSTGterm ratios for n = 80 E-to-I ( g ) and for n = 20 I-to-E ( h ) connections. i – s Six different hundred-neuron networks were simulated with synaptic conductance values drawn randomly from log-normal distributions. Other details were the same as in ( d – h ). i – k Hinton diagrams for one of the networks. l FP and FN rates were computed for all datasets using medianterm filtering without (purple bars) and with deconvolutionterm (green bars). Fractions are out of a population of {Exc. FP: n = 59,520; Exc. FN: n = 480; Inh. FP: n = 59,760; Inh. FN: n = 240} connections. Error bars depict SEM. m Mean f 1 scores for all datasets and connection types. Error bars, SEM. n = 6 networks. p < 0.001, bootstrap test, n = 3500 iterations. n , o eSTGsterm for E-to-I connections plotted against true connectivityterm magnitude (the excitatory conductance, Gie). Colored dots depict n = 80 connections from the dataset shown in ( i ), and gray dots depict n = 400 connections from other datasets. The best linear fit (solid line) and Spearman’s correlation coefficient are shown for the colored dots. p eSTGterm ratios for n = 480 E-to-I connections of all datasets. eSTGterm ratios are negatively correlated with BF. q , r Same as n , o , for I-to-E connections. Colored dots depict n = 20 connections from the dataset shown in ( i ), and gray dots depict n = 100 connections from other datasets. s eSTGterm ratios for n = 120 I-to-E connections of all datasets. Ratios are negatively correlated with BF.

Fig. 6. In real neuronal data, deconvolution-based estimates of spike transmissionterm gain are especially higher in the presence of burstterm spiking activity.

a – c Data were recorded from hippocampal region CA1term of three freely-moving mice using high-density silicon probes. a Example wideband data (band-pass filtered 0.1–7500 Hz) recorded by two sites on the same shank, separated vertically by 60 μm. Orange and brown traces and tick marks correspond to PYR and INT spikes, respectively. b Cumulative distributions of the firing rates of n = 1041 PYR and n = 215 INT (solid lines). Dashed lines show distributions of n = 611 PYR and n = 152 INT used for further analyses ( d , e ). c Distributions of burstterm indices of the same four neuronal groups. d In PYR-to-INT pairs with putative monosynaptic connections, STGterm estimates based on deconvolutionterm are higher than estimates that employ only medianterm filtering. Scatter plot shows eSTGterm dc-median vs. eSTGterm medianterm for all n = 7991 PYR-to-INT pairs in which both eSTGsterm were positive. PYR-to-INT pairs with a consistent monosynaptic CCHterm peak, in which the PYR also exhibits a non-negative burst indexterm·def, are denoted as “analysis pairs” (pink dots; n = 1274 pairs). For the analysis pairs, the eSTGterm dc-median is higher than the eSTGterm medianterm ( p < 0.001, Wilcoxon’s paired signed-rank testterm·stat). e Four example CCHsterm with the medianterm predictor. Insets: The ACHs of the putative presynaptic neurons, which exhibit burstterm spiking. f dcCCHs for the pairs in ( e ). Insets: Conditional rate (CR) histograms based on the CCHsterm (purple) and the dcCCH (green). In all examples, the deconvolution-based eSTGterm is higher than the eSTGterm estimated without deconvolutionterm. g eSTGterm ratios for the full dataset ( n = 1274 pairs), shown as a function of the presynaptic burst indexterm·def. Deconvolution-based estimates are consistently higher than estimates without deconvolutionterm ( p < 0.001, Wilcoxon’s signed-rank test). Estimates are particularly higher when spike bursting is more prevalent ( p < 0.001, permutation testterm·stat).

To focus on the contribution of deconvolutionterm, we chose to use a single timescale separationterm method from this point onward. Of the STG-based methods, medianterm filtering yielded the best performance, with respect to both quantification (Fig. 4h ) and detection (Fig. 4t ). Run times were tested on synthetic spike trains with firing rates of 2 and 8 spk/s, recorded over 5 h. The medianterm [IQR] runtime of medianterm filtering was 0.65 [0.64 0.67] ms, three to four orders of magnitude faster than the runtime of PSP-based methods (GLMCCterm: 3.72 [2.88 4.55] s; CoNNECT: 1.62 [1.54 1.69] s). These runtimes and performance made medianterm filtering our method of choice for evaluating the contribution of deconvolutionterm in biophysical simulationsterm (Fig. 5 ) and in real data (Fig. 6 ).

To investigate how burstterm spiking activity affects STGterm estimation of pairs embedded in a neuronal network and whether deconvolutionterm can improve STGterm estimation, we first conducted fixed-strength simulationsterm. Every E-cell was connected randomly to one I-cell with a fixed excitatory synaptic conductance (Gie = 0.055 mS/cm 2 ), and every I-cell was connected randomly to one E-cell and to one other I-cell with a fixed inhibitory conductance (Gei = 0.25 mS/cm 2 ; Fig. 5d ). Burstiness was randomized. Although E-to-I connection strengths were the same for all pairs, medianterm filtering without deconvolutionterm yielded eSTGsterm which were highly affected by the BF ( ρ = −0.68, p < 0.001, permutation testterm·stat; Fig. 5e ). In contrast, medianterm filtering with deconvolutionterm produced eSTGsterm that were not correlated with burstiness ( ρ = −0.18, p = 0.11, permutation testterm·stat; Fig. 5f ). The ratio between eSTGsterm without and eSTGsterm with deconvolutionterm was close to one when BF was low (−42 mV) and gradually decreased as BF increased ( ρ = 0.891, p < 0.001, permutation testterm·stat, Fig. 5g ). The eSTGterm ratio for the I-to-E connections showed a similar trend, in which the eSTGterm ratio was negatively correlated with the BF ( ρ = −0.835, p < 0.001, permutation testterm·stat; Fig. 5h ). Thus, even when connections have fixed strengths, the eSTGsterm depend on burstterm spiking activity. However, dependency is minimized using deconvolutionterm.

To evaluate the contribution of deconvolutionterm in a more heterogenous setting, we generated networks of conductance-based CA1term model neurons with varied connectivityterm strengths (drawn randomly from log-normal distributions; Fig. 5i ). Other parameters were the same as in Fig. 5d–h . Over six hundred-neuron networks, the mean [SEM] f 1 score obtained with deconvolutionterm was 0.89 [0.008], higher than f 1 without deconvolutionterm (0.84 [0.009]; p < 0.001, bootstrap test; Fig. 5j–m ). In addition to improved detection, deconvolutionterm yielded eSTGterm estimates which were higher than eSTGterm estimates without deconvolutionterm, especially when bursting was prevalent (Fig. 5n–s ). eSTGterm ratios for the E-to-I and I-to-E connections were negatively correlated with BF (E-to-I: ρ = −0.88, p < 0.001, permutation testterm·stat, Fig. 5p ; I-to-E: ρ = −0.837, p < 0.001; Fig. 5s ). Thus, in synthetic neuronal networks, deconvolutionterm improves both the detection and quantification of inter-neuronal connectivityterm.

In real neuronal data, deconvolution-based estimates of spike transmissionterm gain are especially higher when burstterm spiking is prevalent To test the contribution of deconvolutionterm to the analysis of real neuronal data, we obtained data from freely-moving mice. The dataset consisted of 1041 PYR and 215 INT recorded using high-density silicon probes implanted in hippocampal region CA1term of three mice (Fig. 6a ). In these data, the medianterm [IQR] firing rates of PYR and INT were 0.90 [0.14 1.68] spk/s and 9.44 [1.04 17.83] spk/s, respectively (Fig. 6b ). For quantifying burstterm spiking activity of real spike trains, we used the “burst indexterm·def”. Burstterm indices were 0.28 [0.1 0.47] for PYR and −0.283 [−0.37 −0.19] for INT (Fig. 6c ). In contrast to simulationsterm, rSTGsterm are of course unknown for the real neuronal data. Therefore, we tested the contribution of deconvolutionterm by measuring the eSTGterm ratio (without and with deconvolutionterm) as a function of the burst indexterm·def. Only PYR-to-INT pairs with a consistent monosynaptic CCHterm peak ( p < 0.001, Poisson testterm·stat), in which the PYR exhibited a non-negative burst indexterm·def, were analyzed (1274/7991 pairs). For the analysis pairs, firing rates of PYR and INT were 0.96 [0.23 1.7] spk/s and 13.05 [5.27 20.91] spk/s, respectively (Fig. 6b ). Burstterm indices were 0.33 [0.18 0.49] for PYR and −0.28 [−0.348 −0.219] for INT (Fig. 6c ). For the analysis pairs only, deconvolution-based eSTGsterm were consistently higher than eSTGsterm without deconvolutionterm (mediansterm: 0.019 and 0.017, respectively; p < 0.001, Wilcoxon’s paired signed-rank; Fig. 6d ). Thus, in pairs of real neurons with putative monosynaptic connections, deconvolutionterm yields higher eSTGsterm than eSTGsterm estimated without deconvolutionterm.

The higher eSTGsterm yielded by medianterm filtering with deconvolutionterm may result from a fixed DC shift or be burst-dependent. To differentiate between these possibilities, we measured eSTGterm ratios as a function of PYR burst indexterm·def. To visualize the effect of deconvolutionterm, we chose several pairs of units (Fig. 6e, f ) in which the PYR exhibited burstterm spiking activity (Fig. 6e , insets). The CCHsterm between the units exhibited prominent peaks within the causal ROI (0 < τ ≤ 5 ms), consistent with monosynaptic connectivityterm. In addition, the CCHsterm exhibited side lobes on both sides of the main peak, consistent with PYR burststerm (Fig. 6e ). In contrast, the dcCCHs did not exhibit side lobes (Fig. 6f ), and deconvolution-based eSTGsterm were higher than eSTGsterm estimated from the CCHsterm. Over 1274 pairs, eSTGterm ratios exhibited a negative correlation with burstterm indices ( ρ = −0.348; p < 0.001, permutation testterm·stat; Fig. 6g ). Furthermore, most eSTGterm ratios were smaller than one (medianterm: 0.88; p < 0.001, Wilcoxon’s signed-rank test; Fig. 6g ). The negative correlation between eSTGterm ratio and burst indexterm·def in real neuronal data, and the fact that most eSTGterm ratios are smaller than one, are consistent with the controlled simulationsterm (Figs. 3 – 5 ). Together, these observations suggest that deconvolutionterm improves both the quantification and the detection of inter-neuronal connections.

Discussion We showed that the CCHterm between two neurons can be decomposed into a sum of three elements. Based on the mathematical formulation, we developed a deconvolution-based algorithm that removes second-order spike train statistics from the CCHterm, improving connectivityterm quantification and detection. Using pairwise point process and conductance-based network simulationsterm, we showed that burstterm spiking activity modifies the CCHterm, impairing connectivityterm estimates yielded by several distinct methods. Deconvolutionterm removed the effect of burstterm spiking from the CCHterm and recovered the constructed connectivityterm. In noisy point process simulationsterm with known STGsterm, deconvolutionterm reduced the errors of all STG-based methods, and improved detection performance for the tailsterm, medianterm, and GLMCCterm methods. In networks of conductance-based CA1term model neurons, the eSTGsterm yielded by medianterm filtering depended on burstterm spiking activity, and this dependency was removed by deconvolutionterm. In real neuronal data recorded from CA1term region of freely-moving mice, deconvolution-based eSTGsterm were consistently higher than eSTGsterm without deconvolutionterm, in particular when bursting was prevalent. Together, the results imply increased accuracy upon using deconvolutionterm when analyzing real neuronal data.

Estimates of STGterm/PSP yielded by all methods employed (tailsterm, jitterterm, medianterm, GLMCCterm, and CoNNECT) were averaged over time and spike patterns. The underlying assumption was that the ACHs and connectivityterm strengths are fixed over time. However, previous studies have shown that the postsynaptic effect of a presynaptic spike may depend on the presynaptic inter-spike interval 3 , 10 , 11 and on other parameters 32 , 47 . In vitro work showed that synaptic transmission can exhibit short-term plasticity, where the EPSPterm magnitude in a postsynaptic neuron is modified according to the order and timing of presynaptic spikes 48 . In some cases, EPSPterm magnitude decreased over consecutive spikes, and in other cases consecutive spikes caused facilitation of EPSPterm magnitude 3 , 49 . Consistent with the in vitro work, spike transmissionterm estimates derived from the CCHterm in vivo were distinct when CCHsterm were computed using the first, second, or third spike in a presynaptic burstterm 21 , 32 , 50 , 51 . Other in vivo studies used trial-averaging and time-resolved cross-correlation measures, showing that spike-to-spike correlations depend on the time lag relative to a stimulus or an action 32 , 47 , 52 . As is, the deconvolutionterm approach presented here (as well as the timescale separationterm methods employed) does not take into account time-dependent changes in the auto-correlation (e.g., of bursting activity) and cannot estimate time-dependent changes of spike transmissionterm. However, deconvolutionterm together with a timescale separationterm method can be applied to segmented data to track time-dependent changes in spike transmissionterm (as done in refs. 21 , 32 , 47 , 51 , 52 ), yielding estimates differentiated from synchronous firing and burstterm spiking activity.

We showed that deconvolutionterm improves connectivityterm estimates for bursting activity, but deconvolutionterm can remove other effects of second-order spike train statistics on the CCHterm. While bursting is a relatively common deviation from Poisson spiking in real spike trains, other deviations are clearly possible. For instance, in some brain regions, strong periodicityterm is encountered, for instance at the ripple 53 or gamma 54 range. When individual neurons exhibit phase locking to field oscillations, the ACH may exhibit periodicityterm 54 , 55 . Then, the CCHterm between connected neurons will be affected by ACH periodicityterm. Deconvolutionterm can be used to differentiate between the effect of transmitted spikes and the firing patterns of the presynaptic neuron. Thus, deconvolutionterm can minimize the distortion of connectivityterm estimates and improve the analysis of slower CCHterm features.

Following multiple technological advances, the capability to record data from many neuronal pairs has drastically increased. The temporal resolution of calcium imaging techniques has gradually improved 56 – 60 . Furthermore, voltage imaging now allows recording the activity of multiple neurons at sub-millisecond temporal resolutions 61 , 62 . Large-scale electrophysiological techniques have evolved to enable recording hundreds of spike trains from the same or different brain areas simultaneously over prolonged durations 63 – 65 . Analyzing neural circuit dynamics based on massively parallel spike timing data requires adequate quantification of synaptic connectivityterm. In recent years, several analysis methods have been introduced, producing increasingly more accurate circuit maps based on multiple spike train recordings. Deconvolutionterm can improve the outcome of those methods, and has the potential to improve the performance of any CCH-based method. The accuracy and sensitivity improved by deconvolutionterm enables the construction of more accurate neuronal connectivityterm maps which may include weaker or otherwise hidden connections. Detailed maps, which take into account overlooked connectivityterm motifs, may enhance our understanding of the neural mechanisms underlying brain function.

Methods Conductance-based model of a leaky integrate and fire (LIF) neuron driven by a Poisson train To quantify the interplay between spike transmissionterm gain (STGterm), unitary EPSPterm (uEPSP) amplitude, and background activity we constructed a simple conductance-based, leaky integrate and fire (LIF) model neuron with membrane potential V 2 driven by a Poisson spike train via a single excitatory synapse (LIF-syn model). The model was governed by two differential equations and a single condition: C d V 2 t d t = I t − g L V 2 t − E L − g S S ( t ) V 2 t − E S + g N η t 4.1 if V 2 ≥ V th , V 2 ← V reset 4.2 d S t d t = H V 1 1 − S ( t ) / τ r − S ( t ) / τ d 4.3 We used C = 1 μF/cm 2 , g L = 0.5 mS/cm 2 , E L = −60 mV, g S ∈ [0 0.032] mS/cm 2 , E S = 0 mV, g N = 1 mS/cm 2 , V th = −50 mV, V reset = −70 mV, τ r = 0.1 ms, and τ d = 3 ms. Membrane potential variability of the postsynaptic neuron, which may stem from many unknown sources, was quantified in the LIF model by additive Gaussian noiseterm η ( t ) ~ G (0, σ ), with an SD that was varied between runs, σ ∈ [0 5] mV. Whenever a LIF spike occurred, V 2 was held at V peak = 50 mV for T spike = 1 ms before being reset to V reset for another ARP 2 − T spike = 1 ms.

Explicit external inputs to the model were a fixed bias current, I ( t ) = 1 μA/cm 2 ∀ t , and a presynaptic spike train. The presynaptic spike train was generated as a homogeneous Poisson process with λ 1 = 5 spk/s, and then spikes were decimated to maintain a minimal inter-spike interval of 4 ms. Whenever a presynaptic spike occurred, V 1 = V peak = 50 mV for T spike = 1 ms; otherwise V 1 = E L = −60 mV. For incorporating the presynaptic spike train in the synaptic differential equation, we used H ( V ) = (1 + tanh( V /4))/2. Thus, whenever an isolated presynaptic spike occurred, the synaptic variable S rose exponentially (with τ r ) toward one over 1 ms, and then relaxed back to zero (with τ d ). Short-term synaptic plasticity (facilitation and/or depression) was not included. The model was implemented using explicit second-order Runge-Kutta endpoint numerical integration method with a time step of Δ t = 0.1 ms and total duration of T sim = 180 min.

For constructing Fig. 1a , the LIF did not receive any explicit synaptic input ( g S = 0 mS/cm 2 ), but σ ∈ [0 5] mV was varied at 0.5 mV steps (11 levels). For every σ level, six repetitions were performed. For constructing Fig. 1b , we also employed 11 levels of g S , matched to yield uEPSP amplitudes uEPSP ∈ [0 2] mV at 0.2 mV steps. For every { σ,uEPSP } combination, six repetitions were performed.

Quantification of spike transmissionterm gain in the LIF-syn model The firing rate of the LIF, λ 2 , was determined by counting all spikes and dividing by T sim . Denoting the firing rate of the Poisson input by λ 1 , the “real” spike transmissionterm gain (rSTGterm) in the stationary scenario may be intuitively defined as the number of spikes added due to the presence of a synapse, divided by the number of presynaptic spikes, namely △ λ 2 / λ 1 = λ 2 syn − λ 2 no − syn / λ 1 . However, such computation is flawed since the addition of transmitted spikes occasionally decimates spontaneous (noise-induced) spikes. In the case of a LIF model with an absolute refractory period ARP 2 and time step Δ t , the exact rSTGterm is obtained by multiplying the naive ratio by a correction factor that scales the firing rate difference by the fraction of non-occupied bins: rSTGterm = λ 2 syn − λ 2 no − syn λ 1 1 − λ 2 no − syn ⋅ ARP 2 + Δ t . 5 Spike cross-correlation histograms Denote two spike trains by s 1 and s 2 . Each train is a sum of delta functions recorded over a duration T and can be expressed as s 1 t = ∑ i = 1 N 1 δ t − t i , ∀ 0 ≤ t i ≤ T and s 2 t = ∑ i = 1 N 2 δ t − t i , ∀ 0 ≤ t i ≤ T , respectively. Then, the cross-correlation between the two trains at any time lag τ is C C τ = ∫ 0 ∞ s 1 t s 2 t + τ d t . Equivalently, the spike trains can be expressed as ordered lists: s 1 = t i ; i = 1 , … , N 1 , ∀ 0 ≤ t i ≤ T ; and s 2 = r i ; i = 1 , … , N 2 , ∀ 0 ≤ r i ≤ T . Then, the cross-correlation can be expressed as C C τ = ∑ i = 1 N 1 ∑ j = 1 N 2 # t i − r j = τ , where the indicator function is defined as #{ true } = 1 and #{ false } = 0. In practice, two spikes rarely occur at the exact same instant, and the cross-correlation is expressed as a histogram with a finite bin size B and a total of 2 M + 1 bins, spanning the time range from − MB to MB . Then, the CCHterm count at the m th bin ( m ∈ − M , − M + 1 , … 0 , … , M − 1 , M ) sums the cross-correlation over the time range m − 1 / 2 B ≤ τ < m + 1 / 2 B , and is expressed as CCHterm m = ∑ i = 1 N 1 ∑ j = 1 N 2 # m − 1 / 2 B ≤ t i − r j < m + 1 / 2 B . Unless mentioned otherwise, we used a bin size of B = 1 ms and histogram range of 2 M + 1 = 61 ms. We used a bin size of 1 ms to minimize the probability that two spikes of the same train will occur in the same bin, keeping the counting process Poisson. A 1 ms bin size has used in multiple extracellular 22 , intracellular 2 , and across species 66 studies. Furthermore, the CoNNECT method can operate only on CCHterm with bin size of 1 ms.

“Other inputs” predictor In this work we used four methods to estimate the “other inputs” element of the CCHterm: the tailsterm predictor, the jitterterm predictor, the medianterm filter predictor, and the slow part of a generalized linear model for cross-correlations (GLMCCterm). The tailsterm predictor is computed by averaging over all CCHterm bins that are far from the zero-lag bin, both causal and anti-causal bins: pred tailsterm m = 1 K ∑ k CCHterm k , k = − M , − M + 1 , … , − M 0 , M 0 , … , M − 1 , M 6.1 The sum is over all K bins for which | k | ≥ M 0 ; we used M 0 = 11 ms. Therefore, the tailsterm predictor has the same value for all 2 M + 1 CCHterm bins (a straight line). The jitterterm predictor is computed by convolving the CCHterm with a partially hollowed Gaussian kernel: pred jitterterm m = ∑ k = − 3 δ 3 δ CCHterm m − k w k 6.2 The Gaussian kernel, denoted as w [ k ], has a unity sum, standard deviation δ , support of 6 δ , and the central bin is partially hollowed (0.6). Unless noted otherwise, we used δ = 5 ms. Jittering every spike within a rectangular ± δ window centered around the spike (“spike time jitterterm”) many times and computing an average jittered CCHterm converges exactly to the convolutionterm of the CCHterm with a triangular window that has a width of 4 δ + 1 20 . Convolutionterm with a partially hollowed Gaussian window with an SD of δ reduces the false-negative rate and increases detection power 20 . The hollowed medianterm filter predictor is computed by applying a medianterm filter of order 2 δ over the CCHterm, while ignoring the value in the central bin: pred medianterm m = medianterm CCHterm m − δ , … , m − 1 , m + 1 , … m + δ 6.3 The GLMCCterm predictor is given by the slow part, a ( t ), of the GLM defined by refs. 40 , 41 .

Conditional rate CCHterm and spike transmissionterm curve To estimate the spike transmission curve (STC)term·def from the CCHterm, we first compute the conditional rate CCHterm (crCCH). The crCCH is derived by subtracting, for every time lag, the predictor from the CCHterm, and scaling by the number of trigger spikes and the bin size: crCCH m = CCHterm m − pred m / N 1 ⋅ B 7 Thus, the crCCH has units of spikes per second (spk/s). Second, we find the extremum value of the crCCH within the casual temporal region of interest (ROI: 0 < m ≤ 5 ms). The two zero-crossing points, to the left and to the right of the extremum value, are defined as bounds, B L and B R . To preserve causality, if B L falls at the zero or a negative time lag, the first bin after the zero-lag bin is used instead. However, B R may be outside (to the right of) the ROI. The STC is then defined as equal to the crCCH for all bins within the [ B L , B R ] interval, and zero otherwise (Fig. 1c ).

The ROI can be modified according to the expected form of the STC, based on two parameters: the delay and the jitterterm of spike transmissionterm. When the transmission delay is short and the STC is narrow, a short ROI near-zero lag is suitable. Here, an ROI of (0,5] was chosen to ensure including at least part of the monosynaptic STC, consistent with previous results in neocortex 32 and hippocampus exhibiting a transmission delay in the 0–3 ms range 21 , 22 , 67 , 68 .

Computation of spike transmissionterm gain The estimated spike transmissionterm gain (eSTGterm) is the integral over the STC, computed numerically as the sum of all estimated STC values multiplied by the bin size: eSTGterm = ∑ m = B L B R STC m ⋅ B 8 Mathematical framework for deconvolved CCHterm·def The deconvolutionterm approach is based on the realization that the CCHterm between two spike trains can be expressed as an exact sum of four elements and approximated by a sum of three elements. The decomposition (Eq. 17 ) is mathematically exact when the system is assumed to be linear time invariant (LTI). “Linear” implies that the spike transmissionterm process is linear, but does not imply linearity of spike generation, which is a highly nonlinear process. Linearity of spike transmissionterm means that the process of transmitting spikes is additive and does not depend on the history of spike transmissionterm. Linearity is incompatible with neuronal refractoriness, and is assumed for mathematical simplicity. “Time invariant” implies that the impulse responses of spike transmissionterm—and in this specific case, also the ACHs—are constant over the duration of the recording. These assumptions are implicit in the computation of any CCHterm and in most CCHterm applications (see “Discussion”).

Within the LTI framework, we proceed as follows. The specific system studied is a pair of reciprocally-coupled neurons, each of which also receives external input (Fig. 2a ). Formally, the observed spike train of neuron 1, s 1 , is a sum of delta functions and thus has the form s 1 t = ∑ i = 1 N 1 δ t − t i , ∀ 0 ≤ t i ≤ T 9 N 1 is the total number of spikes in s 1 , t i is the instant at which the i th spike occurred, and T is the total recording duration. The spike train of neuron 2, s 2 , is defined in an analogous manner by s 2 t = ∑ i = 1 N 2 δ t − t i , ∀ 0 ≤ t i ≤ T . We denote the finite impulse response of spike transmissionterm from neuron 1 to neuron 2 by h 1 τ , ∀ 0 < τ ≤ τ max 10.1 The impulse response between neuron 2 and neuron 1, h 2 , is denoted in an analogous manner by h 2 τ , ∀ 0 < τ ≤ τ max 10.2 Specifically, h 1 and h 2 are non-zero only for positive time lags, which is the standard definition of causality. Using these definitions, the spike train of neuron 2, s 2 , is the superposition of spikes from two sources: (i) transmitted spikes, due to s 1 filtered by the synaptic connection between the two neurons, h 1 ; and (ii) background spikes, denoted by s 2 0 , which are due to other sources: s 2 t = s 2 0 t + ∫ − ∞ ∞ s 1 t − τ h 1 τ d τ 11 The spike train of neuron 1, s 1 , can be expressed in an analogous manner. To facilitate the analysis, we take the Fourier transform of all quantities, and obtain the frequency domain representation of Eq. ( 11 ): S 2 f = S 2 0 f + S 1 f H 1 f , 12.1 S 2 ( f ) is the Fourier transform of s 2 ( t ), S 2 ( f ) = ∫ − ∞ ∞ s 2 ( t ) exp ( − 2 π i f t ) d t . In an analogous manner, we obtain the Fourier transform of s 1 : S 1 f = S 1 0 f + S 2 f H 2 f 12.2 Equation 12 consists of a pair of coupled equations, for which the solution is: S 1 f = S 1 0 f + S 2 0 f H 2 f 1 − H 1 f H 2 f 13.1 S 2 f = S 2 0 f + S 1 0 f H 1 f 1 − H 1 f H 2 f 13.2 To link the CCHterm with the impulse response, the cross-correlation between the two spike trains CCHterm 12 τ = ∫ − ∞ ∞ s 1 t s 2 t + τ d t 14.1 is written in the frequency domain as F CCHterm 12 = S 1 f ¯ S 2 f 14.2 where S 1 f ¯ denotes complex conjugation, due the fact that when computing convolutionsterm—but not cross-correlations—one of the signals is time reversed. Plugging Eq. (13) into Eq. ( 14.2 ) we obtain F CCHterm 12 = S 1 0 f ¯ S 2 0 f + S 1 0 f ¯ S 1 0 f H 1 f + S 2 0 f ¯ S 2 0 f H 2 f ¯ + S 2 0 f ¯ S 1 0 f H 2 f ¯ H 1 f 1 − H 1 f H 2 f 2 15 The elements in the numerator of Eq. ( 15 ) can be rewritten in terms of correlation functions of the latent (background) spike trains as follows: CCHterm 12 0 τ = S 1 0 f ¯ S 2 0 f = ∫ − ∞ ∞ s 1 0 t s 2 0 t + τ d t 16.1 CCHterm 21 0 τ = S 2 0 f ¯ S 1 0 f = ∫ − ∞ ∞ s 2 0 t s 1 0 t + τ d t 16.2 ACH 1 0 τ = S 1 0 f ¯ S 1 0 f = ∫ − ∞ ∞ s 1 0 t s 1 0 t + τ d t 16.3 ACH 2 0 τ = S 2 0 f ¯ S 2 0 f = ∫ − ∞ ∞ s 2 0 t s 2 0 t + τ d t 16.4 Yielding F CCHterm 12 = F CCHterm 12 0 + F ACH 1 0 H 1 f + F ACH 2 0 H 2 f ¯ + F CCHterm 21 0 H 2 f ¯ H 1 f 1 − H 1 f H 2 f 2 17 This completes the exact decomposition of the CCHterm into four additive elements in the frequency domain. When the overall gain of the feedback loop is small, i.e., when H 1 ( f ) H 2 ( f ) → 0, the last element and the denominator of Eq. ( 17 ) vanish, and the Fourier transform of the CCHterm simplifies to F CCHterm 12 = F CCHterm 12 0 + F ACH 1 0 H 1 f + F ACH 2 0 H 2 f ¯ 18.1 Which is represented in the time domain as CCHterm 12 = s 1 0 ⋆ s 2 0 + ACH 1 0 * h 1 + ACH 2 0 * h 2 − τ 18.2 In Eq. ( 18.2 ), ⋆ denotes cross-correlation and * denotes convolutionterm. We have reached an important conclusion. When the overall gain of the feedback loop is small, the cross-correlation between two spike trains can be expressed as a sum of three elements: (1) the cross-correlation between the uncoupled spike trains; (2) the convolutionterm of the auto-correlation of one uncoupled spike train with the impulse response of the coupling to the second train; and (3) the convolutionterm of the auto-correlation of the second uncoupled train with time-reversed coupling to the first train.

The exact same derivation applies to inhibitory connections, the difference being that the plus signs in Eq. ( 11 ) and Eq. (12) are replaced by minus signs. In the absence of connectivityterm, Eq. ( 17 ) simplifies to CCHterm = CCHterm 12 0 . Alternatively, in the presence of purely feedforward connectivityterm, e.g., when H 2 f ¯ → 0 , Eq. ( 17 ) simplifies to F CCHterm 12 = F CCHterm 12 0 + F ACH 1 0 H 1 f 19.1 Which is represented in the time domain as CCHterm 12 = s 1 0 ⋆ s 2 0 + ACH 1 0 * h 1 19.2 In some systems, e.g., when neurons are strongly coupled 69 , 70 , the approximation H 1 ( f ) H 2 ( f ) → 0 cannot be made. In general, the approximation H 1 ( f ) H 2 ( f ) → 0 can be made when the loop gain (rSTGterm 12 times rSTGterm 21 ) is lower than a minimal unidirectional STGterm, which can be detected with finite data. The detection of a connection depends on multiple variables including recording duration, firing rates, and the method employed. Therefore, an exact number for a minimal detectable unidirectional STGterm cannot be obtained. However, in the specific case of unidirectional connection with an impulse response spanning a single bin and Poisson firing pattern of both neurons, the minimal detectable unidirectional STGterm can be computed exactly using the tailsterm predictor. Under these conditions, the CCHterm baseline is 20 λ = F 1 F 2 TB , where F 1 and F 2 are the firing rates of the pre- and postsynaptic neurons in spk/s, T is the effective recording duration, and B is the CCHterm bin width (both measured in seconds) For a given λ and alpha (detection threshold) level, the inverse Poisson distribution can be used to obtain the minimal number of transmitted spikes required for detection. The minimal STGterm is then calculated as the minimal number of transmitted spikes, divided by F 1 T . For example, the minimal STGterm which can be detected for F 1 = 1 spk/s, F 2 = 10 spk/s, at T = 50, 000 s using B = 0.001 s , and an alpha level of 0.001 is 0.00142.

Algorithm for deconvolving ACHs from the CCHterm In general, we do not have access to the latent variables s 1 0 and s 2 0 . In other words, we cannot determine which spikes in a train, e.g., s 2 , are background spikes (e.g., s 2 0 ), and which are transmitted from s 1 (Eq. 11 ). Thus, in general, ACH 1 0 and ACH 2 0 are unknown and Eq. 18 cannot be evaluated. However, if the number of transmitted spikes is small, compared to background spikes, then the ACHs of the background trains can be estimated from the measured spike trains: ACH 1 0 τ ≈ ACH 1 τ = ∫ 0 ∞ s 1 t s 1 t + τ d t 20.1 ACH 2 0 τ ≈ ACH 2 τ = ∫ 0 ∞ s 2 t s 2 t + τ d t 20.2 Where the lower integration bound is 0 rather than −∞ due to the non-negativity of the time axis. The CCHterm is estimated directly from the spike trains as in Eq. ( 14.1 ), with the time axis modified in the same manner: CCHterm 12 τ = ∫ 0 ∞ s 1 t s 2 t + τ d t 20.3 Finally, we define H ~ 1 f = H 1 f F ACH 2 21.1 H ¯ ~ 2 f = H 2 f ¯ F ACH 1 21.2 Using these measurements (Eq. 20) and definitions (Eq. 21) in Eq. ( 18.1 ) and rearranging, we obtain H ~ 1 f + H ¯ ~ 2 f + F CCHterm 12 0 F ACH 1 F ACH 2 = F CCHterm 12 F ACH 1 F ACH 2 22 Where the r.h.s. is derived from the observed spike trains. Taking the inverse Fourier transform and defining the “other inputs” element as OI = F − 1 F CCHterm 12 0 F ACH 1 F ACH 2 23 We obtain the basis of the deconvolutionterm algorithm: h ~ 1 + h ~ 2 − τ + OI = F − 1 F CCHterm 12 F ACH 1 F ACH 2 24.1 Due to causality (Eq. 10), h 1 and h 2 (− τ ) reside on opposite sides of the zero lag. The r.h.s. of Eq. ( 24.1 ) is denoted the “deconvolved CCHterm·def” (dcCCH): dcCCH 12 = F − 1 F CCHterm 12 F ACH 1 F ACH 2 24.2 When unidirectional deconvolutionterm is applied, Eqs. ( 24.1 – 2 ) can be replaced by: h 1 + OI = F − 1 F CCHterm 12 F ACH 1 24.3 dcCCH 12 = F − 1 F CCHterm 12 F ACH 1 24.4 Once the dcCCH is determined numerically (see below), the “other inputs” element can be determined using timescale separationterm (a predictor; Eq. 6). Following subtraction and rescaling (Eq. 7 ), the STCs and the eSTGsterm (Eq. 8 ) are evaluated.

The first step in the actual deconvolutionterm algorithm is to compute the two count ACHs and the CCHterm from the spike trains. The second step is to scale the ACHs, such that each of the ACHs will have a sum of one. To scale the count ACH, the zero-lag bin is first set to zero. Then, the mean is subtracted, every bin in the zero-mean ACH is divided by the total number of spikes, and the zero-lag bin is set to complement the sum to one. Scaling ensures that deconvolutionterm will be affected only by the shape of the ACHs, and not by the number of spikes in each train. If the count ACH is exactly flat, the normalization yields a delta-like histogram (one at zero-lag and zero everywhere else), and deconvolutionterm does not modify the CCHterm. Third, the r.h.s. of Eq. (24) is evaluated.

The deconvolutionterm algorithm as derived above (Eq. 20–24) is based on Eq. 18, but makes two key simplifying assumptions. First, that the shape of the ACH is not considerably distorted by the transmitted spikes (Eq. 20). Second, the impulse response is not considerably distorted by the divisive scaling factor inherent in the derivation (Eq. 21).

With respect to Eq. 20, using point process simulationsterm, we found that even when the presynaptic neuron spikes at a relatively high rate (10 spk/s) with strong bursting activity (BR = 0.4), and the postsynaptic neuron spikes at a lower rate (3 spk/s), the estimates of connectivityterm are not distorted as long as the rSTGterm is below 0.25. Notably, strong connections and high firing rates are not characteristic of local cortical networks. However, different settings may appear in different systems, e.g., thalamocortical feedforward connections 28 . When connectivityterm is very strong and presynaptic firing rates are high, prior knowledge about connectivityterm can be used to employ unidirectional deconvolutionterm (Eq. 24.3 –4) instead of the “standard” deconvolutionterm (Eqs. 24.1 – 2 ).

While ACHs effects are eliminated by the deconvolutionterm algorithm, a divisive scaling factor (the denominator of Eq. 21) is introduced which can distort the recovered impulse responses. If connectivityterm is unidirectional and the postsynaptic train is Poisson, the scaled ACH 2 is a delta function, the divisive factor in Eq. ( 21.1 ) is unity, and h ~ 1 = h 1 . Furthermore, if connectivityterm is unidirectional, then even if the presynaptic train is non-Poisson, the numerator of the r.h.s. of Eq. ( 21.2 ) is zero and thus h 2 = 0. When connectivityterm is bidirectional and one of the trains is non-Poisson, or when connectivityterm is unidirectional and the postsynaptic train is non-Poisson, the divisive factor may become consequential. Then, prior knowledge about connectivityterm can be used to employ unidirectional deconvolutionterm (Eq. 24.3 –4).

Point process spike train simulationsterm For investigating different STGterm estimation methods, we devised a point process simulationterm that generates two spike trains of a two-neuron network with precisely defined STGsterm. Each of the synthetic spike trains, s 1 ( t ) and s 2 ( t ), was initially generated by random sampling from a distinct rate function, λ 1 ( t ) and λ 2 ( t ), respectively, and then spikes were added or removed according to predetermined STCs. Thus, simulationsterm consisted of three sequential steps: (1) generation of a rate function for each neuron; (2) stochastic sampling of spike trains, independently from each rate function; (3) modification of each spike train according to the other spike train and the corresponding STC. All simulationsterm were generated over N = T /Δ t samples, where Δ t = 0.001 s is the time step and T is the duration of the simulationterm.

To generate a rate function for one neuron, we first set a desired mean firing rate λ d which was then modified based on the desired gamma order γ and burst fractionterm·def BR according to λ m = λ d ∙ γ /(1 + BR ). For instance, for generating a spike train with gamma order γ = 2 the firing rate is doubled, and for a spike train in which half of the spikes are first spikes in a two-spike burstterm the firing rate is halved. For generating a fixed rate function, we set λ ( t ) = λ m ∀ t ∈ [1 N ], separately for each neuron. To create co-modulated spike trains, we generated correlated rate functions for the two neurons. For that, the rate function of each neuron was time dependent, defined as: λ ( t ) = λ m + λ m ∙ λ c ( t ). The co-modulation rate function λ c is a clipped pink noiseterm signal, generated by filtering Gaussian white noiseterm sampled at Δ t intervals, η ( t ) ~ G (0, σ c ), with a decaying exponential e − 1 / τ c , where τ c = 20 ms. The pink noiseterm was then clipped to the [−1 1] range, yielding the zero-mean λ c ( t ). The same co-modulation function was used for generating the two rate functions, λ 1 ( t ) and λ 2 ( t ).

In the second step, each spike train was generated by stochastic sampling from the corresponding rate function. For each sample, we determined whether a spike did or did not occur, according to random sampling from a Binomial distribution with parameters B ( n , p ) = (1, λ ( t )Δ t ). When λ ( t ) is time-varying or constant, the resulting spike trains correspond to non-homogenous or homogenous Poisson processes, respectively. After generating a spike train, three modifications were made. First, spikes were removed according to the gamma order γ : to realize an n th-order gamma process, every n th spike was retained, and all other spikes were decimated. Retaining all spikes corresponds to a first-order gamma (i.e., Poisson) process. Second, burstterm spiking activity was added by inserting spikes based on the burst fractionterm·def ( BR ) parameter. For a non-zero BR , each spike was defined to be a “first in burstterm” randomly with probability BR 1 . The precise timing of the “second in burstterm” spike was determined by random (uniform) sampling from a 5-element symmetric triangular window that lagged the first spike in the burstterm by two samples. A “third in burstterm” spike was added with probability BR 2 ; then, the overall BR is BR 1 + BR 1 ∙ BR 2 , and the timing of the third spike was drawn from a three-sample lagged window. In actual simulationsterm, BR 1 was varied in the [0 0.4] range, whereas BR 2 = 0.4 was held constant. Third, ARP considerations were imposed by removing all spikes that occurred shortly (ARP = 2 ms) after another spike.

In the third step, connectivityterm was added between the two spike trains. The same STC waveform was used in all simulationsterm: a 5-element asymmetric triangular window that lagged the presynaptic spike by one sample. However, the STC integral (the STGterm) was varied in sign and value, mimicking excitatory and inhibitory connectivityterm of various gains. For each spike of the presynaptic neuron, the occurrence and timing of the added (or decimated, for inhibitory connections) spike was determined by random (uniform) sampling from the STC. For excitatory transmission, this procedure can generate more than one postsynaptic spike, as required since the STGterm is not limited to the [0 1] range. Connectivityterm was implemented independently in every direction. After applying connectivityterm between the two trains, ARP considerations were imposed once more.

For creating noisy point process simulationsterm (Fig. 4 ), we simulated 1250 spike train pairs. In 60% of the pairs, co-modulation dynamics were added. For the simulationsterm with co-modulation, the co-modulation parameter σ c was drawn randomly from a uniform distribution between 1 and 15 spk/s, with a fixed co-modulation time constant ( τ c = 20 ms). The firing rates of the pre- and postsynaptic neurons were 2 and 8 spk/s, respectively. The burst fractionterm·def of the presynaptic neuron was varied uniformly between 0 and 0.4. Recording duration was set randomly between 90 min and 5 h. Of the 1250 simulated pairs, 500 pairs were simulated with monosynaptic excitatory connections; 500 pairs were simulated with monosynaptic inhibitory connections; and the other 250 pairs were simulated without connectivityterm. Connectivityterm strengths (rSTGterm) were drawn randomly from a log-normal distributionterm·stat 43 – 45 with {mean, SD} of {0.019, 0.01} for the excitatory connections, and {−0.014, 0.007} for the inhibitory connections.

Conductance-based model of CA1term network To model synaptic connectivityterm within a network of neurons with ACHs and CCHsterm resembling those observed in real neuronal data (Fig. 6 ), we generated a network of conductance-based E- and I-cells. The E-cell model, derived from the “simple model” of Izhikevich 71 , included dynamics on the membrane potential ( V ), and a slower, phenomenological recovery variable ( u ). In addition, the model included synaptic input and noiseterm. The model equations are: C d V d t = I in e t + V k V − E L e V − V t − u − I synaptic + g N η t d u d t = U a U b V − E L e − u if V > V peak then V ← V reset ; u = u + U step 25 Membrane potential variability, which may stem from many unknown sources (e.g., unrecorded neurons), was modeled by an additive noiseterm term, generated by random sampling from a zero-mean Gaussian distribution η ( t ) ~ N (0, σ ) independently for every cell. In this model, spike waveforms were modeled explicitly by the quadratic in Eq. 25 . Whenever a spike occurred (i.e., V peak was crossed), the membrane potential V was reset to V reset , and the recovery variable u was incremented by U step . By modifying these two parameters, a given E-cell model was tuned to fire predominantly single spikes, spike burststerm, or mixes thereof (Fig. 5a ). The specific parameters values used for the E-cell model are detailed in Table 1 .

Table 1.

Parameters used for modeling E-cells in the conductance-based networks.

Parameter Value Units C 1 μF/cm 2 V k 0.01 mS/(cm 2 ·mV) E L e −60 mV V t −45 mV g N 1 mS/cm 2 U a 0.02 mS/cm 2 U b 0.01 mS/cm 2 V peak 40 mV V reset [−42 −35] mS/cm 2 U step [1.12 1.225] mV·mS/cm 2 σ e 2 mV I in e 0 μA/cm 2 N e For the I-cells, we used the Wang-Buzsáki model 72 , describing the dynamics of the membrane potential ( V ), sodium inactivation ( h ), and delayed-rectifier potassium ( n ). The full model also included synaptic currents and noiseterm, and reads C d V d t = I i n i t − g L V − E L i − g N a h m ∞ V 3 V − E N a − g K n 4 V − E K − I synaptic + η t + g N η t d h d t = h ∞ V − h τ h V d n d t = n ∞ V − n τ n V 26.1 The gating variables for the I-cell ( x = h,m,n ) had voltage-dependent time constants ( τ x ) and steady-state values ( x ∞ ) as follows: h ∞ V = 0.07 e − V + 58 20 0.07 e − V + 58 20 + 1 1 + e − V + 28 10 , τ h V = 0.2 0.07 e − V + 58 20 + 1 1 + e − V + 28 10 26.2 m ∞ V = 0.2 V + 35 1 − e − V + 35 10 0.2 V + 35 1 − e − V + 35 10 + 4 e − V + 60 18 26.3 n ∞ V = 0.01 V + 34 1 − e − V + 34 10 0.01 V + 34 1 − e − V + 34 10 + 0.125 e − V + 44 80 , τ n V = 0.2 0.01 V + 34 1 − e − V + 34 10 + 0.125 e − V + 44 80 26.4 Other parameters values used for the I-cell model are detailed in Table 2 . Synaptic connections were modeled as in ref. 73 . For the e ’th E-cell, the total synaptic current was I synaptic , e = ∑ j = 1 N e g e j S e j V e − E s e , j + ∑ k = 1 N i g e k S e k V e − E s i , k 27.1 Where N e ( N i ) is the number of E-cells (I-cells). The notation g ej indicates the maximal synaptic conductance from presynaptic E-cell j to postsynaptic E-cell e . All excitatory-to-excitatory (E-to-E) synapses had the same reversal potential, regardless of the presynaptic neuron ( E se,j = E se , ∀ j ), but the maximal excitatory conductance value g ej could vary between presynaptic neurons. In practice, we set all g ej = g ee = 0. All inhibitory-to-excitatory (I-to-E) synapses had the same reversal potential, regardless of the presynaptic neuron ( E si,k = E si , ∀ k ). The maximal inhibitory conductance g ek was either fixed for all I-cells ( g ek = g ei ; Fig. 5d ) or varied between neurons (Fig. 5i ). All synaptic activation variables corresponding to the same presynaptic neuron had the same dynamics, regardless of the postsynaptic neuron ( S ej = S j , S ek = S k , ∀ e ). For the i ’th I-cell, the total synaptic current was modeled by I synaptic , i = ∑ j = 1 N e g i j S i j V i − E s e , j + ∑ k = 1 N i g i k S i k V i − E s i , k 27.2 Table 2.

Parameters used for modeling I-cells in the conductance-based networks.

Parameter Value Units C 1 μF/cm 2 g L i 0.1 mS/cm 2 E L i −65 mV g Na i 35 mS/cm 2 E Na i 55 mV g K i 9 mS/cm 2 E K i −90 mV g N 1 mS/cm 2 σ i 2 mV I in i −0.75 μA/cm 2 N i All excitatory-to-inhibitory (E-to-I) synapses had the same reversal potential ( E se,j = E se , ∀ j ), but the maximal conductance values, g ij , could be varied between E-cells. All inhibitory-to-inhibitory (I-to-I) synapses had the same reversal potential ( E si,k = E si , ∀ k ), but the maximal conductance values g ik could be varied between I-cells. All synaptic activation variables corresponding to the same presynaptic neuron had the same dynamics ( S ij = S j , S ik = S k , ∀ i ).

For an excitatory/inhibitory presynaptic neuron, the dynamics of the corresponding synaptic variable ( S e / S i ) depended on the presynaptic membrane potential ( V e / V i ) and the synaptic rise and decay time constants, following: d S e d t = H V e 1 − S e τ r e − S e τ d e 28.1 d S i d t = H V i 1 − S i τ r i − S i τ d i 28.2 H V = 1 + tanh V / 4 / 2 28.3 All synaptic parameters values used are detailed in Table 3 . Numerical integration was done using the explicit second-order Runge-Kutta endpoint (modified Euler) method with integration time step of Δ t = 0.025 ms and simulationterm duration of T sim = 240 min.

Table 3.

Synaptic parameters used in simulating conductance-based networks.

Parameter Value Units Notes τ r e 0.1 ms AMPA τ d e 3 ms AMPA E e 0 mV AMPA τ r i 0.3 ms GABA A τ d i 9 ms GABA A E i −80 mV GABA A g ie 0.055 mS/cm 2 Fig. 5i : mean: 0.055, SD: 0.02 g ee 0 mS/cm 2 g ei 0.25 mS/cm 2 Fig. 5i : mean: 0.22, SD: 0.08 g ii 0.25 mS/cm 2 Fig. 5i : mean: 0.22, SD: 0.04 Animals and ethics Three freely-moving male C57BL/6J mice were used in this study (Table 4 ). At the time of implantation, mice aged 14–16 weeks and weighed 26–30 g. After implantation, mice were single-housed to prevent damage to the implanted apparatus. All animal handling procedures were in accordance with Directive 2010/63/EU of the European Parliament, complied with Israeli Animal Welfare Law (1994), and approved by the Tel Aviv University Institutional Animal Care and Use Committee (IACUC #01-16-051).

Table 4.

List of experimental animals.

Animal ID Sex Age a [week] Weight a [g] Probe Sessions All PYR b All INT b Analysis PYR c Analysis INT c Analysis pairs c Rank CC d P -value e mDS1 Male 14 25.7 Dual-sided64 9 322 50 170 37 427 −0.19 0.0003 mP30 Male 14 28.6 Stark64 8 277 57 184 39 387 −0.3 0.0003 mP31 Male 16 30 Stark64 10 442 108 257 76 458 −0.09 0.0776 Summary 27 1041 215 611 152 1272 −0.35 0.0003 a At the time of implantation.

b “All PYR” and “All INT” refer to the total numbers of well-isolated units recorded during the listed sessions.

c “Analysis PYR”, “Analysis INT”, and “Analysis pairs” refer to the subsets used for actual analyses (cf. Fig. 6b–d ).

d “Rank CC” refers to Spearman’s rank correction coefficient between eSTGterm ratios and burstterm indices.

e “ P -value” refers to permutation testterm·stat comparing the rank CC to a zero null.

Probes and surgery Each animal was implanted with a multi-shank silicon probe (Diagnostic Biochips) attached to a movable micro-drive. The probes used were Stark64 (two mice) and Dual-sided64 (one mouse). The Stark64 probe consists of six shanks, spaced horizontally 200 µm apart, with each shank consisting of 10–11 recording sites, spaced vertically 15 µm apart. The Dual-sided64 probe consists of two dual-sided shanks, spaced horizontally 250 µm apart, with each shank consisting of 16 channels on each side (front and back), spaced vertically 20 µm apart. In all mice, probes were implanted in the neocortex above the hippocampus (PA/LM, 1.6/1.1 mm) under isoflurane (1%) anesthesia 38 , 74 . After every recording session, the probe was translated vertically downward by up to 70 µm. The recorded data included only recordings from the CA1term pyramidal cell layer, recognized by the appearance of multiple high-amplitude spiking units and iso-potential spontaneous ripple events.

Recording procedures Neuronal activity was recorded in 4.4-h sessions (medianterm of 27 sessions; range, 3–12 h). Every session started with a baseline neural recording of at least 15 min, while the animal was in the home cage or in a 0.8 m diameter open field. After the baseline recordings, the animal ran on a linear track, received optogenetic stimuli, or both, for a period of at least 60 min. Sessions ended with another baseline period of at least 30 min. Only data recorded during spontaneous activity, in the lack of any optogenetic stimuli, were used in this work.

Spike detection and sorting Neural activity was filtered, amplified, multiplexed, digitized on the headstage (0.1–7500 Hz, x192; 16 bits, 20 kHz; RHD2132 or RHD2164, Intan Technologies), and recorded by an RHD2000 evaluation board (Intan Technologies). Offline, spikes were detected and sorted into single units automatically using either KlustaKwik3 75 , 76 or KiloSort2 77 . Automatic spike sorting was followed by manual adjustment of the clusters. Only well-isolated units were used for further analyses (amplitude >40 µV; L-ratio <0.05 78 ; ISI index <0.2 79 ). Units were classified into putative pyramidal cells (PYR) or parvalbumin-immunoreactive [PV]-like inhibitory interneurons (INT) using a Gaussian mixture modelterm·stat 29 .

Quantifying burstterm spiking activity For quantifying burstterm spiking activity of real spike trains, we used a “burst indexterm·def” 80 . The precise definition employed was burst indexterm·def = head − tailterm head + tailterm , where head is the sum of all ACH counts in the 2 < τ ≤ 10 ms range, and tailterm is the sum of all ACH counts in the 35 < τ ≤ 50 ms range.

Selection of a subset of data For evaluating the contribution of deconvolutionterm to the analysis of real neuronal data, we used spike trains of units recorded from CA1term of three freely-moving mice during spontaneous activity. We denote PYR-to-INT pairs for which both eSTGterm medianterm and eSTGterm dc-median were positive as the “All pairs” group (7991/9030 simultaneously-recorded PYR-INT pairs; Fig. 6d , blue dots). For eSTGterm comparisons, we focused on a subset of pairs (the “Analysis pairs” group; 1274/7991 pairs; Fig. 6d , pink dots) composed of pairs which met four criteria. (1) The trigger unit (PYR) exhibited a non-negative burst indexterm·def (Fig. 6c ). (2) An excitatory connection was detected by the medianterm filter with and without deconvolutionterm. (3) The ratio between eSTGterm medianterm and eSTGterm dc-median was between 0.5 and 2 (28 pairs were excluded due to this criterion). (4) The total number of counts in the count CCHterm (at the −30 < τ ≤ 30 range) was above 400 (5 pairs were excluded due to this criterion).

Detection analysis For quantifying detection performance of each method, we first computed false-positive (FP) and false-negative (FN) rates for each type of connection. For excitatory connections, the FP rate was defined as the number of times in which the method detected an excitatory connection when no connection was simulated, divided by the total number of tested connections ( N conn ). The excitatory FN rate was defined as the number of times in which an excitatory connection was simulated but was not detected by the method, divided by N conn . For inhibitory connections, the FP rate was defined as the number of times in which the method detected an inhibitory connection when no connection was simulated, divided by N conn . The inhibitory FN rate was defined as the number of times in which an inhibitory connection was simulated but was not detected by the method, divided by N conn . Next, we computed the f 1 score, defined as: f 1 = TP TP + 0.5 ( FP + FN ) 29 Where TP is the true positive rate, which is defined as the number of times in which a simulated connection (excitatory or inhibitory) was correctly detected by the method, divided by N conn . The FP and FN used in Eq. ( 29 ) are the sum of the excitation and inhibition FP and FN rates: FP = FP exc + FP inh and FN = FN exc + FN inh , respectively.

To determine whether the f 1 score of one group (group 1 ) is consistently larger than the f 1 score of another group (group 2 ), we devised a resampling with replacement (bootstrap) testterm·stat. For each resampled dataset, we first computed the f 1 score for each group, and then compared the two f 1 scores. We repeated the process 3500 times. We then counted the fraction of times for which group 1 f 1 score was smaller than group 2 f 1 score. Based on the resampled scores, we computed the mean and SEM.

Slope comparison To determine whether the linear slope of one group (group 1 ) is consistently larger than the slope of another group (group 2 ) we devised another, more conservative bootstrap test. For each resampled dataset, we fitted a slope for each group. We repeated the process 300 times, resulting in 300 different slopes for each group and 90,000 slope pairs. We then counted the fraction of pairs for which group 1 slope was smaller than group 2 slope. Pairing slopes from different resampling iterations provides a more conservative test than only pairing slopes from the same iteration.

Statistics and reproducibility In all statistical tests used in this study, a significance threshold of αterm·stat = 0.05 was used. Two exceptions were the threshold used for determining whether two units exhibit monosynaptic connectivityterm, and the classifier used to determine whether a unit is a PYR or an INT ( α = 0.001). For Fig. 5m , reproducibility was carried out across networks (Fig. 5m ); for Fig. 6g , reproducibility was carried out across animals (Table 4 ). All descriptive statistics ( n , mean, medianterm, range, IQR, SEM) can be found in the results and figure legends. To examine whether a group medianterm is larger or smaller than an expected value, we used Wilcoxon’s signed-rank test (one-tailed)term·stat. Differences between two group mediansterm were tested with Mann–Whitney’s U -test for unpaired samples (two-tailed)term·stat, or with Wilcoxon’s signed-rank paired test for paired samples (two-tailed). Association between parameters was estimated with Spearman’s rank correlation and tested using a permutation testterm·stat.

Reporting summary Further information on research design is available in the Nature Research Reporting Summary linked to this article.

Supplementary information Peer Review File (4.7MB, pdf) 42003_2022_3450_MOESM2_ESM.pdf (471KB, pdf) Description of Additional Supplementary Files Supplementary Data 1 (13.6MB, xls) Reporting Summary (1.2MB, pdf) Acknowledgements This work was supported by ERC #679253, ISF #638/16, Rosetrees #A1576, and CIHR-IDRC-ISF #2558/18. We thank Ortal Amber-Vitos, Refaela Atsmon, Roni Gattengo, Tom Manovitz, Shir Sivroni, and Alexander Tarnavsky-Eitan for their constructive comments.

Author contributions L.S., A.L., S.S., and E.S. built optoelectronic devices; L.S. and E.S. implanted animals, developed analyses, and wrote the paper; L.S., A.L., and H.S. collected data; L.S. analyzed data and prepared figures; E.S. developed the mathematical framework and the deconvolutionterm algorithm, conceived, and supervised the project.

Peer review Peer review information Communications Biology thanks Soumya Jana, Anuththara Rupasinghe and the other, anonymous, reviewer(s) for their contribution to the peer review of this work. Primary Handling Editors: Enzo Tagliazucchi and Gene Chong. Peer reviewer reports are available.

Data availability The source data behind the graphs in the paper are available in Supplementary Data 1 . The data used in this study are available from the corresponding author upon reasonable request.

Code availability Code for implementing deconvolutionterm and medianterm filtering, along with example spike data, are hosted publicly on GitHub, accessible via https://github.com/EranStarkLab/CCH-deconvolution .

Competing interests The authors declare no competing interests.

Footnotes Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information The online version contains supplementary material available at 10.1038/s42003-022-03450-5.