Inferring monosynaptic connections from paired dendritic spine Ca2+ imaging and large-scale recording of extracellular spiking

Objective: Techniques to identify monosynaptic connections between neurons have been vital for neuroscience research, facilitating important advancements concerning network topology, synaptic plasticity, and synaptic integration, among others. Approach: Here, we introduce a novel approach to identify and monitor monosynaptic connections using high-resolution dendritic spine Ca2+ imaging combined with simultaneous large-scale recording of extracellular electrical activity by means of high-density microelectrode arrays. Main results: We introduce an easily adoptable analysis pipeline that associates the imaged spine with its presynaptic unit and test it on in vitro recordings. The method is further validated and optimized by simulating synaptically-evoked spine Ca2+ transients based on measured spike trains in order to obtain simulated ground-truth connections. Significance: The proposed approach offers unique advantages as (a) it can be used to identify monosynaptic connections with an accurate localization of the synapse within the dendritic tree, (b) it provides precise information of presynaptic spiking, and (c) postsynaptic spine Ca2+ signals and, finally, (d) the non-invasive nature of the proposed method allows for long-term measurements. The analysis toolkit together with the rich data sets that were acquired are made publicly available for further exploration by the research community.


Introduction
Monosynaptic connections enable the flow of information between neurons. The contact point at which transmission occurs is the synapse, where the postsynaptic component of excitatory inputs is often a small protrusion called 'dendritic spine' . The extent to which individual excitatory synapses contribute to the drive that the postsynaptic neuron receives depends on a multitude of synaptic properties, such as synaptic strength (with pre-and postsynaptic components), short-term plasticity processes, and the synapse location within the dendritic tree [1][2][3]. Moreover, the local dendritic environment, including adjacent synaptic inputs and dendritic properties, may allow for synaptic interactions and local dendritic information processing [4,5]. The precise presynaptic spike patterns are, in this context, often of crucial importance. In addition, presynaptic activity and its relationship to postsynaptic spiking shapes the development of the synapse over time [6]. An in-depth interrogation of the functional properties, of organizational principles, and of the development of dendrites and their inputs requires, therefore, advanced methodologies, which provide both structural and functional information, ideally over extended periods of time.
Existing approaches for the identification and probing of monosynaptic connections have specific advantages and drawbacks [7,8]. For example, axonal tracing in fixed tissue may allow for the exact localization of the synapse but cannot be used for functional investigations, while dual patch-clamp recordings grant, for a short period of time, access to pre-and postsynaptic activity, but the location of synaptic contacts is not directly available. Another approach to infer functional connectivity relies on correlated activity of extracellularly recorded spike trains [9][10][11]. Here, the pairwise cross-correlation of unit spike trains is computed, based on which the spike-transmission probability (STP) between preand postsynaptic neuron can be assessed. While the non-invasive nature of this approach is certainly a strong advantage, STP values for a given connection depend on synaptic integration with other inputs and, consequently, they depend on the network state [12]. This implies that some connections may not be identifiable based on spike trains alone. Additional drawbacks of the spike-time cross-correlation approach include that an identification of excitatoryto-excitatory connections often requires particularly long recordings of spiking activity [13], and that information about postsynaptic subthreshold activity is not available. Finally, methods that provide exquisite control over synaptic activation, either by electrical stimulation or using optical means (e.g. glutamate uncaging or optogenetics), do not provide access to precise spontaneous presynaptic spike times.
Here, we introduce a multimodal approach to map monosynaptic connections in functional neuronal networks that fills some of the existing gaps. To do so, we exploit the high temporal resolution of extracellular recordings and the high spatial resolution of Ca 2+ imaging. Using simultaneous highdensity microelectrode array (HD-MEA) recordings and subcellular-resolution Ca 2+ imaging in rat primary neuronal cultures, we were able to record precise spiking sequences of the many potential presynaptic cells and the single-spine Ca 2+ transients of the evoked postsynaptic signals. We developed a novel and robust analysis pipeline to map monosynaptic connections on the basis of these paired recordings. We validate our method using data-driven simulations and demonstrate with several paired experimental recordings that the method can indeed be used to detect monosynaptic connections.
The presented approach represents a versatile tool for functional characterization of synaptic connections, including synaptic plasticity, activitydependent synapse development and synaptic integration. In addition, location-dependencies of various synaptic properties can be probed.
All the data sets acquired in this work are available in the Neurodata Without Borders format [14] on the DANDI archive https://gui.dandiarchive.org/ #/dandiset/000223. The source code to implement the method and notebooks to reproduce the figures are available at: https://github.com/starquakes/mea_ spineca_mapping. A complementary-metal-oxide-semiconductor (CMOS) based HD-MEA featuring 26 400 electrodes was used [15]. This HD-MEA system has 1024 reconfigurable readout channels that can be simultaneously recorded from at 20 kHz sampling frequency. We packaged the HD-MEA chip with epoxy (Epo-Tek 353ND, 353ND-T, Epoxy Technology Inc. Billerica, MA, United States) to cover the gold bond wires and the printed circuit boards before use. We used large plastic rings (diameter of 35 mm) for maintaining the neuronal cultures on the chips so as to enable access and sufficient scanning movement of the high-magnification dipping lenses (40×, 60×) of the upright confocal microscope. The electrodes at a pitch of 17.5 µm were covered with electro-deposited platinum black to decrease electrode impedance and improve the signal-to-noise characteristics.

HD-MEA system and electrode selection
The electrode configurations for simultaneous readout were selected in two ways: (a) a brief prescan of all electrodes was performed, and the active electrodes were subsequently chosen. When there were more than 1024 active electrodes, we occasionally generated two to three sets of configurations. (b) Alternatively, we first identified the position of the target cell on the chip and then selected the surrounding electrodes.

Animal use and cell culture preparation
All experimental protocols involving animals were approved by the Basel-Stadt veterinary office according to Swiss federal laws on animal welfare, and were carried out in accordance with the approved guidelines. Prior to culturing cells, the HD-MEA was sterilized for 40 min in 70% ethanol and rinsed three times with de-ionized (DI) water. Next, the electrode array was treated with 20 µl of 0.05% (v/v) poly(ethyleneimine) (Sigma-Aldrich) in borate buffer (Thermo Fisher Scientific, Waltham, MA, United States) at 8.5 pH, for 40 min at room temperature. This coating improves cell adhesion and makes the surface more hydrophilic. The chips were then rinsed three times with DI water to remove the remaining solution. Next, we added 8 µl of 0.02 mg ml −1 laminin (Sigma-Aldrich) in Neurobasal medium (Gibco, Thermo Fisher Scientific) to support the growth and differentiation of the cells. The chips with laminin were incubated for 30 min at 37 • C. Cortices of E-18 Wistar rat embryos were dissociated in trypsin with 0.25% ethylenediaminetetraacetic acid (EDTA) (Gibco) and counted with a hemocytometer as described previously [16]. We then seeded 15 000-20 000 cells on top of the electrode arrays. Subsequently, the chips were incubated at 37 • C for 30 min before adding 2 ml of plating medium. The plating medium stock solution consisted of 450 ml Neurobasal (Invitrogen, Carlsbad, CA, United States), 50 ml horse serum (HyClone, Thermo Fisher Scientific), 1.25 ml Glutamax (Invitrogen), and 10 ml B-27 (Invitrogen). After three days, 50% of the plating medium were replaced by a growth medium, with stock solutions consisting of 450 ml Dulbecco's Modified Eagle Medium (D-MEM) (Invitrogen), 50 ml Horse Serum (HyClone), 1.25 ml Glutamax (Invitrogen), and 5 ml sodium pyruvate (Invitrogen). The procedure was repeated twice a week. The chips were kept inside an incubator at 37 • C and 5% CO 2 . All the experiments were conducted between days in vitro (DIV) 18 and 21.

Fluorescence imaging of dendritic spine Ca 2+ signals
For confocal fluorescence microscopy, we used a Nikon NiE upright microscope, equipped with Yokogawa W1 spinning disk scan head, an Andor iXon Ultra EMCCD camera (Oxford Instruments), and 40×/0.80 NA or 60×/1.00 NA waterobjectives (Nikon). Cells were transduced with floxed jGCaMP7b (AAV1-syn-FLEX-jGCaMP7b-WPRE; Addgene #104493, MOI = 5 × 10 5 vg) [17] and lowtiter Cre (AAV9-hSyn-Cre-WPRE-hGH; Addgene #105553, MOI = 5 × 10 3 vg) adeno-associated viruses (AAVs) on DIV 7. We used a sparse expression of the calcium indicator to reduce the risk that the measured spine Ca 2+ trace would be contaminated by the activity of other neurites. The use of confocal microscopy further reduced the risk of such contaminations. We expect that our approach would still work with a non-sparse calcium indicator expression, under the condition that some neurons exhibit isolated dendritic branches without adjacent neurites from other cells. A sparse expression, however, is likely to maximize the number of spines that can be tested for a connection. A customized clamping device had parafilm M attached to the objective and the culture to avoid water evaporation during long recording sessions. Confocal image time series containing the target postsynaptic cell were acquired at 7-15 Hz, and large-scale electrical recordings of extracellular signals were simultaneously performed with the HD-MEA setup. A typical recording length was 2-3 min. Synchronization was achieved by sending transistor-transistor logic (TTL) pulses from the camera to the field-programmable gate array controller of the HD-MEA system for post-recording alignment.

Processing of Ca 2+ fluorescence traces
Fluorescence traces were initially extracted by averaging over the pixel values in the region of interest (ROI) for each frame. Next, the ∆F/F was computed by extracting the baseline value with a moving percentile function (25th percentile, 100 data points as the moving stride). This also effectively detrended the signal. In order to isolate synaptically-evoked Ca 2+ transients from the fluorescence traces, we removed the contamination by Ca 2+ signals associated with backpropagating action potentials (bAPs) of the imaged cell. To achieve this, a robust linear regression was performed using the dendritic shaft and spine traces and the contribution of the dendritic shaft (i.e. the bAP-component) was subtracted from the spine trace [18].

Processing of extracellular recordings
The extracellular recordings were first spikesorted using SpyKING Circus 0.8.4 (width of templates = 3 ms, spike threshold = 6, cut-off frequencies for band-pass Butterworth filter = 300 Hz/9500 Hz, spatial radius considered = 210 µm), followed by manual curation of each recording. Templates with a similar extracellular signature were deemed likely to be recorded from the same neuron and the spike time cross-correlograms were employed to support merging.

Burst detection based on global firing rate
To detect periods of strong synchronized network bursting, we analyzed variations in the population firing rate. As a first step, we computed individual unit spike rates. Next, hyperactive units were excluded to ease the computation of burst boundaries later in the procedure. To do so, we followed two strategies: (a) to exclude units whose rates were beyond the 99.75th percentile of the distribution of unit firing rates or (b) to exclude units whose rates were above 15 Hz.
The population rate of the clean units was computed in 10 ms bins. The p th baseline and p th thresh percentiles of these rates were assigned as the baseline and threshold rates, respectively ( f baseline and f thresh ). Peaks in the population rate trace that satisfied the following conditions were then identified as bursts: peak rate was at least f thresh and its prominence was at least f thresh . Prominence, a measure of how much a peak stands out relative to its neighbors, was computed as the height of the peak relative to its lowest local valley [19]. Additionally, a minimal inter-peak interval of t ipi was imposed.
Individual burst boundaries were computed in the neighborhood of each peak using a smoothed version of the network rate trace. Smoothing was performed using locally weighted linear regression over a span of n span data points. We found that smoothing yielded more robust estimates of burst boundaries. The start and end of each bursting period was defined as the first time point at which the smoothed rate crossed f baseline to the left and right of the peak.
Edge cases, where data began or ended during a burst, were excluded. Further, a minimal inter-burst interval of t ibi was imposed. If bursting periods were closer than t ibi , they were merged into a single larger burst.
The following parameter values were used in our analysis: p baseline = 85, t ipi = 50 ms, n span = 5, t ibi = 75 ms. Depending on the degree of synchrony observed in the data, p thresh was visually chosen for each data set. Lower values resulted in a less stringent definition of a bursting period, while higher values imposed more stringent conditions, resulting in much smaller periods being marked as bursts. Our choice of p thresh typically lays between the 96th and the 99th percentile. Finally, we extended the detected bursting periods by 0.1 s before the burst onset and 0.5 s after the burst, in order to account for the slow dynamics of the Ca 2+ signal.

Simulation of synaptically-evoked spine Ca 2+ transients
We used a set of release probabilities following a lognormal distribution [20] (figure 2(b)). We used the numpy.random.lognormal function to generate a set of log-normally distributed release probabilities with the mean value of 2.5 and the standard deviation of 0.6 and scaled it by 0.01 to shift the values towards the desired range. Some very few values above 1 were replaced by random values between 0 and 1. We selected a background ROI close to the spine to measure the noise values and used the spine ROI as the signal region. In order to estimate the signal-to-noise ratio (SNR), the same number of pixels were extracted from the noise region and the signal region. The standard deviations of individual pixel traces were calculated, and we used the median values of the standard deviations of the signal and noise regions as the amplitudes of the simulated signal and noise, respectively. Then, the simulated spine response could be written as: where C is the baseline of the fluorescence trace, λ denotes the amplitude of the pure signal (difference between the standard deviation of the measured signal and recording noise), k(t) is the instantaneous firing rate of the spike train, and (k * φ)(t) denotes the temporal convolution of k(t) with the exponential kernel φ(t) (decay τ = 0.5 s). σε(t) indicates the recording noise-σ is the amplitude of the recording noise and ε(t) the white noise [21].

Binary classification
As described in the results section, we performed binary classification using the outcomes from correlation tests with and without the ground-truth units. A set of quantile thresholds ranged from the 0th to 100th quantile (figure 4(b)) of the two distributions of GT-Rs and test-Rs. Note that in rare cases, the best-R of the correlation test with ground-truth unit was not from the ground-truth unit. For such an outcome, we classified it as false positive (FP) if it passed the thresholds, or true negative (TN) if it did not. Then the F-score was calculated as follows: Here we used 0.5 for β, as precision weighs more than recall in the test. The optimal thresholds could be identified when the F-score reached its maximum.

Identifying the presynaptic unit of a dendritic target spine
The proposed method for associating a specific dendritic spine with its putative presynaptic unit relies on a 'brute-force' approach (figure 1). We record the spiking activity of a large fraction of a neuronal network by means of HD-MEAs, while we simultaneously image the synaptically-evoked spine Ca 2+ signals. Subsequently, we compute for each of the recorded spike-sorted units an approximation of the spine Ca 2+ trace that would be evoked in a postsynaptic neuron of the respective unit. Finally, a pairwise Pearson correlation between the measured spine Ca 2+ transients and the computed traces is used to identify the most likely presynaptic unit. For data acquisition, an in-house-developed CMOS HD-MEA recording unit [15] with 26400 electrodes at 17.5 µm pitch, featuring up to 1024 channels that can be simultaneously recorded, was mounted under a spinning-disc upright confocal microscope ( figure 1(a)), as the HD-MEAs are not transparent. Primary rat cortical neurons were plated exclusively onto the active sensing area of the chip, and sparse expression of the fluorescent calcium indicator was achieved by co-transduction of floxed jGCaMP7b [17] and Cre with two AAVs. Experiments were performed at DIV 18-21. The HD-MEA system allows for selecting an arbitrary configuration of the 1024 electrodes that can be simultaneously used for recording, and different selection strategies are possible. Typically, our selection was based on a pre-scan of all 26 400 electrodes in order to determine electrodes that exhibit spiking activity (see section 2 for details).
In figure 1(b), a more detailed overview of the processing pipeline is shown. Following temporal alignment, the HD-MEA and Ca 2+ imaging data are further preprocessed in two separate streams. For Ca 2+ imaging data, the first step requires a selection of the ROI containing the target spine (figure 1(c)), followed by selection of an ROI containing the adjacent dendritic shaft (additional examples of Ca 2+ indicator expression are shown in figure 1-figure supplement 1). The extracted Ca 2+ traces are subsequently ∆F/F transformed. At this stage, the spine trace is a composition of Ca 2+ signals associated with synaptic activations, but also bAPs. Since our method relies on correlations of synaptically-evoked signals, bAP-associated Ca 2+ is an undesired component. Therefore, it is subsequently removed in a demixing step using the bAP-Ca 2+ -dominated shaft trace, as described before [18]. An example of the original ∆F/F of a spine and its adjacent dendritic shaft trace, as well as the demixed spine trace after removing the bAP-Ca 2+ component, is shown in figure 1figure supplement 2. The preprocessing of HD-MEA data includes two main steps: spike sorting and convolution. Spike sorting was performed using SpyK-ING Circus 0.8.4 [22] and SpikeInterface [23], followed by manual curation (using SpyKING Circus Matlab GUI). Next, we employed a decaying exponential kernel [24] that fitted the kinetics of the Ca 2+ indicator [17] (decay τ = 0.5 s) to convolve each of the sorted spike trains (from now referred to as convolved traces). Finally, we downsampled the convolved traces to match the frame rate of the Ca 2+ imaging. Examples of three convolved and downsampled traces are shown in figure 1(d).
Since we aim to use correlations to infer synaptic connectivity, strong bursts of synchronous action potentials in the network will affect the possibility of identifying the presynaptic unit for a given spine. Moreover, increased non-linear postsynaptic interactions [25] during network bursting periods could introduce additional variation between measured and convolved spine Ca 2+ traces. Therefore, an accurate detection and exclusion of network bursts is important for our analysis. By exploiting the network spiking information, we implemented a burst detection algorithm based on a global firing rate principle (see section 2) and excluded periods of strong bursting from the analysis (for an example, see figure 1-figure supplement 3).
As the final step to identify the presynaptic unit for the target spine, we computed the Pearson's correlation between the demixed spine trace and the convolved traces from all recorded units (excluding the bursting periods). The recorded units could then be sorted by their Pearson correlation coefficient (Pearson's R) ( figure 1(e)), and the best matching one (from now referred to as best-matching unit) was identified as the candidate to form a monosynaptic connection to the target spine. Note the remarkable correlation of presynaptic spike times and the convolved trace (blue) for the best-matching unit and the demixed measured spine Ca 2+ trace (orange)especially for smaller events. Moreover, the R value for the best-matching unit (best-R) shows the largest increase with respect to the next (ordered) R value among all units.

Acceptance test for putative monosynaptic connections
In the previous section, we outlined the experimental approach and analysis pipeline to identify a presynaptic candidate unit. However, (a) an incomplete coverage of the network or (b) synchronous firing of neurons could result in a unit displaying the strongest correlation, although it is not the presynaptic neuron of the target spine. The latter is a possibility because presynaptic release probabilities and recording noise will cause stochastic deviations between measured demixed and convolved trace. Therefore, a more formal procedure to accept only a 'trustworthy' presynaptic candidate unit and its putative connection is required. Here we introduce a novel data-driven approach, based on the surrogate method, to perform such an acceptance test.
We propose that each putative connection needs to meet two criteria to be accepted. These two criteria assure that only the spike train of a specific candidate presynaptic unit can, with reasonable confidence, be assumed to generate the measured spine Ca 2+ transients. First, we require the best-R value to fall within a distribution of possible R values for the best-matching unit, assuming biologically plausible presynaptic release probabilities (P r ) and realistic recording noise. To build such a distribution, we construct simulated ground-truth connections of the best-matching unit, i.e. we generate surrogate spine Ca 2+ traces evoked by the best-matching unit spike train (here referred to as ground-truth surrogates; figure 2). As an initial acceptance test, the 95th quantile of the test-Rs (dashed magenta line) and 10th quantile of the GT-Rs (dashed blue line) were selected to determine threshold values that the best-R value (black line) was required to exceed for a putative monosynaptic connection to be accepted.
Essentially, this distribution can be obtained by, again, convolving the unit spike train with the exponential decay kernel that mimics the spine Ca 2+ response. However, this time, we introduce two layers of uncertainty. The first one is given by the release probability of the synapse ( figure 2(a)). Here, we assume a log-normal distribution of P r values that is consistent with previous reports [20]. For each surrogate generation for a given 'presynaptic' spike train, a P r value is drawn from this distribution and then determines the probability of the spikes that are retained as simulated successful synaptic transmissions. Several noise sources during Ca 2+ microscopy introduce additional variations in the recorded spine traces. To account for these variations, we also estimate the SNR for each target spine ( figure 2(b)) and add accordingly scaled white-noise to our surrogate traces (see section 2 for details). Following these steps, an arbitrary number of spine Ca 2+ trace surrogates can be generated for a given presynaptic spike train ( figure 2(c)). For our first acceptance criterion, we generated 1000 spine trace surrogates based on the spike train of the best-matching unit and then calculated for each surrogate the Pearson's R with the convolved trace of the best-matching unit (blue ground-truth surrogate data in figure 3). In this way, we generated a distribution of possible Rs (GT-Rs) for our candidate presynaptic unit. We required the best-R value to fall well within the GT-R distribution, as will be further detailed below.
This first criterion was introduced to assure that Ca 2+ traces that could be evoked by the bestmatching unit were compatible with our experimentally observed best-R value. However, to be able to trust this putative connection, we next needed to exclude that the best-R value also would be a likely outcome for any of the other units. This next step could be reduced to conducting a test for the unit with the second largest R (second-best matching unit) of the correlation between the demixed spine Ca 2+ trace and convolved traces. Essentially, we asked here: are the data compatible with the case that the secondbest matching unit would be connected to our target spine? For this test, we again generated a population of 1000 surrogate spine responses, but this time, we used the spike train of the second-best matching unit (magenta test surrogate data in figure 3). Also for these surrogates, we subsequently calculated the distribution of Pearson's R (test-Rs) by correlating each test surrogate with the convolved trace of the best-matching unit. Our second criterion for the acceptance of a putative monosynaptic connection was then that the best-R did not fall well within the test-Rs distribution. This means that the spike train of the second-best matching unit was unlikely to evoke spine Ca 2+ responses that were compatible with our experimentally observed best-R value. Examples of test-R and GT-R distributions, together with the best-R values for this recording, are shown in figure 3(b). With these distributions, we could now perform the acceptance test based on our two criteria. Finally, it remains to be formalized what is meant by the best-R falling (not) well within a distribution. Empirically, one could choose different quantiles in order to obtain the threshold values for the GT-R and test-R value distributions (for example, in figure 3(b) we used the 95th quantile for test-Rs and the 10th quantile for GT-Rs). In the next section, we will show how to determine data-driven optimized quantile thresholds that maximize the method's performance.

Method validation and threshold optimization
As part of the acceptance test, described in the previous section, we applied fixed quantiles (95th/10th) to the two R distributions (test-R/GT-R) in order to obtain thresholds for the best-R value. However, for best performance, these quantiles should be tailored to each network, where factors such as network synchrony may influence the optimal quantile thresholds. The acceptance test can be considered a classification task, and the best threshold values are those that achieve a compromise between precision and recall. To obtain the optimal quantile thresholds for the acceptance test and to characterize the best performance that can be achieved by our method, we further extended the surrogate method described in the previous section. For each recorded unit i, we generated N surrogate spine Ca 2+ response traces. Each of these presynaptic units could then be considered to be part of a simulated ground-truth connection, and the surrogates constituted simulated spine Ca 2+ measurements ( figure 4(a)). Next, we computed the GT-Rs and test-Rs distributions and we performed the acceptance test with a grid search over quantile values for the test-Rs (Q1) and the GT-Rs (Q2) ( figure 4(b)). From the outcome of the test over all surrogates and all units, we could compute the number of true positives (TPs), when unit i was accepted as a monosynaptic connection, and false negatives (FNs), when unit i was rejected as a monosynaptic connection (or another unit was found). In addition, we ran the test after removing the presynaptic unit i, based on which the surrogates were generated, from the list of recorded units. In this case, the test could yield TNs, when the putative connection was rejected, or FP, when a monosynaptic connection was accepted despite the fact that the presynaptic unit i was  4(a)). After counting all TPs, FNs, TNs, and FPs, we computed the F-score for each pair of Q1-Q2 as a performance metric. Figure 4(c) shows the heat map of F-scores for one recording, with the green arrow showing the maximum score and, therefore, the optimal thresholds for this specific network.

Performance evaluation on experimental data
To evaluate the performance of the proposed method upon application to real experimental data, we conducted 19 paired spine Ca 2+ imaging and HD-MEA recording experiments in ten different cell cultures (table 1 shows the full analysis results). First, the optimal threshold values for the acceptance test were computed for each recording using large-scale simulations, based on the spike-sorted unit data, as described in the previous section. The F-scores associated with the optimal threshold values provided a characterization of the optimal performance. We obtained typically large F-scores across experiments (mean = 0.692 ± 0.136 based on 19 recordings with a mean of 102 ± 43 units per recording and 1000 generated surrogate spine Ca 2+ response traces per unit), which supported our approach and implied that the network characteristics (e.g. synchrony) permitted a good classification in most cultures. Next, we applied our method to measured spine Ca 2+ traces. From each recording, we then typically selected one medium-large dendritic spine that exhibited signs of activity (in one of recordings two spines were chosen). For 9 out of 20 spines, the best-matching unit passed the acceptance test, and a monosynaptic connection was confirmed. An acceptance rate of approximately every other tested connection indicated that a sufficient number of connections for statistical analyses could be obtained a reasonable time period.
We continued our investigation by probing in more detail the basis underlying the acceptance or rejection of a connection. Although the accepted connections had significantly larger best-R values (0.584 ± 0.149) compared to the rejected ones (0.331 ± 0.124) (Mann-Whitney U = 12.0, n 1 = 9, n 2 = 11, P < 0.05), it is important to note that the value of the best-R alone was not a predictor for the outcome of the test. In fact, as one can see from figure 5(b), the acceptance test rejected connections with R values as high as ∼0.6, while in another case an R value of ∼0.2 was enough for acceptance. Figure 5(e) shows two examples of accepted connections, while figure 5(f) displays two rejected ones. Note that for some recordings the test threshold was larger than the GT threshold and vice versa for other recordings. Indeed, we found that putative connections were rejected in three different scenarios: the best-R was smaller than either one of the two thresholds or both thresholds (see summary of all tested spines in figure 5(a)). This finding demonstrates the performance advantage of having two thresholds in the acceptance test. Finally, we further examined the hypothesis that network activity characteristics would affect the performance of our method. We calculated two complementary measures for each recording. First, we determined the ratio of spike events occurring outside and inside of network burst periods (OUT/IN column in table 1), which captured to what extent neuron spiking was entrained by network bursting. Secondly, we calculated the spike contrast of the network, a fast measure of spike train synchrony [26] (implemented in the Elephant  figure 5(c)) of a network correlated with the classification performance as measured by the maximum F-score of the respective network. We have, therefore, ways to screen for the most promising cultures, in terms of identifying monosynaptic connections, before even beginning with spine Ca 2+ imaging.

Discussion
In this paper, we introduced a novel approach for mapping monosynaptic connections using subcellular resolution Ca 2+ imaging of dendritic spines and simultaneously recorded extracellular electrical spiking activity from large-scale networks by means of HD-MEAs. We established a data-driven analysis pipeline that is able to match an optically recorded spine trace with the spiking activity of the putative presynaptic unit and accepts or rejects the connection based on a surrogate strategy. We validated and assessed the performance of the proposed method using simulated ground-truth monosynaptic connections and showed how the developed method can provide a robust detection of monosynaptic connections in experimental data. Our evaluation was performed with 2D-cell-culture recordings. Indeed, the possibility to control the number of cells in the network and the fact that HD-MEA systems can simultaneously record from large fractions of the grown networks, renders dissociated 2D cultures, plated on the HD-MEA chips, particularly suitable for the developed connection mapping method. Recent advances in cell reprogramming technologies have opened up new ways to model neurological diseases by means of induced pluripotent stem cellderived neurons obtained from patients. Many of the modeled diseases are associated with impaired synaptic connectivity [28] and experiments, as detailed below, conducted in patient-derived cultures, could shed light on some of the earliest synaptic dysregulations. In addition, newly developed CMOS HD-MEA systems for in vivo applications [29][30][31] hold the promise that the developed method will also be useful in a 3D context, and a first evaluation could include the application of our optimization and validation pipeline to in vivo HD-MEA data. This way, one could assess the effect of in vivo network characteristics (e.g. synchrony) on the connection mapping performance.
Next, we consider some limitations of the proposed approach and try to discuss possible improvements. In our experimental in vitro data, we achieved an approximately 50% success rate in finding and accepting the presynaptic unit of a target spine. The main potential reasons for rejection are that the presynaptic cell either (i) was not captured by the selection of recording electrodes, (ii) showed little inter-network burst activity, or (iii) exhibited strongly correlated activity with other neurons. Point (i) can easily be addressed by performing multiple paired recordings from the same target spine, while sequentially changing the set of extracellular recording electrodes, so that neurons across the entire HD-MEA surface can be probed (e.g. see figure 1(a)). Alternatively, next-generation CMOS HD-MEAs, featuring an increased number of simultaneously usable read-out channels (∼20 000) [32,33], will provide an even more comprehensive coverage of the network. The other two reasons for rejection are also important to consider. We selected target synapses based on the spines showing at least some individual activity in a pre-recording period, while there was also a significant number of spines that either remained silent or were active mostly during network bursting periods. Therefore, addressing points (ii) and (iii) may be particularly beneficial in order to maximize the number of identifiable connections. The required desynchronization of network activity could be achieved, for example, by applying temporarily desynchronizing stimuli or pharmacological treatments or by using different development stages of the cell cultures. Moreover, the activated states of in vivo neural activity may be naturally conducive to the identification of connections [34]. Finally, some of the assumptions underlying our method could be adjusted to more precisely fit different experimental conditions. For example, we here assume a certain log-normal distribution of release probabilities as reported before [20]. However, the user may wish to incorporate extra information about a specific network. Our analysis pipeline, therefore, allows for a convenient adjustment of the release probability distribution. Moreover, we assumed that spine Ca 2+ transients can be approximated by convolution of the presynaptic spike train with a simple exponential decay kernel without a rising phase. Such simple assumptions have been shown to perform very well for spike deconvolution methods based on Ca 2+ imaging data [24]. Nonetheless, the performance of our connection mapping method may be improved by explicitly modeling the effect of potential synaptic interactions and short-term plasticity on the evoked Ca 2+ transients. However, as we excluded periods of strong network bursting, where these phenomena are most likely to occur, the extra computational costs may not necessarily justify a presumably small gain in performance. In the current implementation, we used Ca 2+ indicators to readout postsynaptic spine signals. Alternatively, the same approach could be applied using next-generation voltage indicators [35] to provide an even higher temporal resolution, which might make the correlations between spike trains and spine signals more robust (provided that an SNR comparable to Ca 2+ can be achieved). Another possible improvement to increase the throughput of identified connections concerns the selection of target spines. Currently, spines that are tested for a connection need to be manually selected. A fully automated selection is a challenging task, but a future integration of such algorithms will drastically reduce the analysis time.
We envision that the proposed method will facilitate various studies of dendrite organization, function, and development ( figure 6). Here we highlight some of the possible investigations that exploit the unique set of provided information, including precise spike patterns of the presynaptic units and of large parts of the network, in conjunction with a postsynaptic Ca 2+ readout. The spatial organization of synapses is critical for the input-output transformation of a neuron [1][2][3]. It is now possible to address questions such as: are there any location dependencies concerning features of presynaptic spiking patterns (e.g. firing rate, bursting behavior, and more)? Furthermore, following the identification of the postsynaptic unit, it is possible to calculate the STP between pre-and postsynaptic neurons. How are STP values organized across the dendritic tree during spontaneous network activity? Such investigations are not possible with existing methodologies. STP estimates, for example, are extracted from spike train cross-correlations that require extremely precise spike times-inference of spike times from Ca 2+ imaging data would not suffice in this context [36][37][38]. Exploiting the combined information of presynaptic spike times and evoked postsynaptic events provided by our method, neurotransmission and its modulation by presynaptic short-term plasticity [2] processes could now be examined in detail at individual synapses during spontaneous presynaptic spike patterns. Typically, a measured synaptic Ca 2+ response trace is associated with transmission failures while the convolved trace represents the response trace with a release probability equaling one. Appropriate scaling of the convolved trace can be achieved by matching of Ca 2+ responses to individual presynaptic action potentials. A direct comparison of the measured Ca 2+ trace and the scaled convolved trace then reveals how the probability of transmission failures changes with time, and it is thus possible to assess if presynaptic shortterm plasticity occurred during certain events (e.g. comparing periods of high vs. low network activity). Moreover, combined long-term HD-MEA recordings with additional periodic spine imaging can be used to address important questions concerning activity-dependent synapse development, by assessing, e.g. structural spine plasticity or changes in neurotransmission. Another possible research direction could include the probing of synaptic interactions. Those are conventionally investigated by glutamate Top: activity-dependent synaptic development. By combining long-term HD-MEA recordings with periodic spine imaging, it is possible to assess if synaptic features (e.g. the size of the target spine (orange) or presynaptic properties (blue)) change with, for example, presynaptic activity levels. Depicted is a putative homeostatic down-regulation of synaptic strength with increasing presynaptic spike rate. Bottom-right: spatial organization of synaptic attributes. One can now examine whether there is a (synapse-soma) location-dependence of properties, such as STP (extracted from unit spike train cross-correlations). Even with only few identified connections for individual postsynaptic neurons, the pooling of data from multiple cells will enable to identify general patterns (e.g. differences between proximal and distal synapses). Bottom-left: short-term synaptic plasticity. The deviation of the (demixed) measured spine Ca 2+ trace from the convolved trace of the best-matching unit reveals how presynaptic release probability varies with activity levels. In the depicted example, heightened levels of presynaptic spiking appear to be associated with an increased probability of transmission failures compared to low-activity periods, which is consistent with short-term synaptic depression.
uncaging experiments [25,39], where presynaptic spiking information is not available. Our method could potentially enable researchers to identify multiple adjacent synaptic inputs and, subsequently, allow to examine interactions during spontaneous activity. In summary, we developed a versatile experimental setup and analysis pipeline for the identification and interrogation of monosynaptic connections that complements existing approaches.

Data availability statement
All the figures and results can be reproduced using the notebooks available on: https://github.com/star quakes/mea_spineca_mapping.