1.

Introduction

Spatiotemporal dynamic patterns of local neural circuits are essential for discovering the mechanisms of the nervous system.^{1, 2, 3, 4} Recently, with the advances in calcium indictors and the two-photon microscopy technique, reconstructions of the action potential (AP) discharge in local neural circuits from their simultaneously recorded calcium fluorescence transients in vivo and in vitro were carried out.^{5, 6, 7, 8, 9} Furthermore, two-photon fluorescence random scanning microscopy has been designed and constructed, through which calcium fluorescence transients sampled at about a 1-kHz frequency from neuronal local circuits can be recorded synchronously.^{10, 11, 12} These data at a high sampling rate enable us to reconstruct AP firing of local neuronal circuits at a millisecond-scale precision.¹³ Although the neuronal burst discharge carries more coding information than a single AP firing,^{14, 15, 16} revealing neuronal AP burst trains from the corresponding calcium waves is well-known as hard work.^{17, 18, 19, 20} This is so not only because burst firing is caused by the large amount of depolarization^{15, 16} and has not been clearly defined, but also because neuronal calcium fluorescence transients always have a nonlinear relationship with their APs firing, causing the diversities and complexities of calcium fluorescence waves.²¹ Practically,^{14, 15, 16} some burst firing sequences consist of APs with intervals of less than 20 ms.

Driven by the necessity to reveal the principles of neuronal circuits, many valuable methods to reconstruct AP sequences from their corresponding calcium fluorescence transients have been designed. For example, the template-matching method was developed to detect a single AP discharge from a calcium fluorescence transient, which works well whenever neurons discharge sparsely.^{22, 23} The uncorrected²⁴ and corrected first-derivative method including deconvolution^{7, 25, 26} can provide reliable neuronal firing rates, but are limited to inferring neuronal AP trains from calcium fluorescence transients with a low signal-to-noise ratio (SNR), defined as the ratio of the amplitude of the template to the standard deviation of noise in the initial calcium waves. Some statistical methods such as clustering²⁷ and a combination of principal component analysis and a support vector machine,²⁸ can only partially describe the calcium wave diversities that derive from the neuronal irregularly discharging at a high frequency rate. Some sophisticated methods^{9, 29, 30} are suitable for reconstructing the AP train from a burst firing. However, compared to the methods just mentioned, more setting parameters or hypotheses may confine their applications. Recently, a method termed “peeling” was proposed for this purpose, but a fixed template matching the rising duration is still necessary.¹³

For the calcium fluorescence transients at a high frequency (1 kHz), the nonlinear calcium dynamics triggered by neuronal burst firing have become increasingly obvious and are difficult to describe using a single model or several patterns. In addition, the SNR of fluorescence signals is a decisive factor for the accuracy of the reconstruction. Therefore, aimed to reconstruct the AP burst firing trains from calcium fluorescence transient with relative low SNRs, a method based on multiscale filter combined with template convolution (MSF_TC) is proposed. First, the parts of calcium fluorescence transients corresponding to neuronal AP firing are identified using a multiscale filter, and the AP train candidates are generated. Then, we combine the AP train candidates with the template function, depending on the fitting degrees with the original calcium waves, to confirm AP firing timing sequences.

To illustrate the MSF_TC algorithm performance, this method was applied in processing real data sets from our optical system.^{11, 12} The reconstruction results of the proposed method and that of two typical methods, the template-matching method and the deconvolution method, were compared. The reconstruction results were further analyzed to emphasize some important features of the MSF_TC algorithm.

2.

Methods

2.2.

Reconstructing Neuronal AP Train

The MSF_TC algorithm contains three steps. The first step is the extraction of the template, which is then used in the following two steps. The AP train candidates are generated in the second step. Last, through template convolution, the optimal AP train is confirmed and regarded as the reconstruction results.

2.2.1.

Extraction of the template

The template extraction is primarily based on the fact that the calcium wave triggered by single AP remains almost unchanged.^{32, 33} Thus, 10 to 20 calcium waves from individual calcium fluorescence transients are aligned at the time when they are ascending most rapidly, and then averaged. The averaged trace is then regarded as the template. For demonstration in Figs. 1 to 1 a template is used to be convoluted with the AP train. Note that the value of the template baseline should be rectified to zero.

Fig. 1

MSF_TC algorithm for AP train reconstruction. (a) calcium fluorescence transients after 10- and 100-Hz low-pass filtering, respectively. (b) Identify the valid calcium fluorescence transient and extract its rising duration [the box in (b)]. (c) The extracted rising duration (middle) corresponds to 5 APs (top), and the predicted APs firing at the peak of the differential of the rising duration of the calcium waves, as indicated by the arrows (bottom). Note that the peak should be more than the given threshold (horizontal line in the bottom panel). Five predicted APs can generate 32 AP train candidates according to Eq. 4. (d) and (f) Two of 32 AP train candidates that convolved with the template and simulated signals were generated. The AP train is regarded as a reconstructed AP train, which corresponds to the simulated signal that has the most similar real signal. Based on this, the AP train in (d) is the reconstructed AP train rather than that in (f).

2.2.2.

AP train candidates generation

In this step, identification of a valid calcium fluorescence transient and the extraction of its rising duration are necessary. The 10- and 100-Hz low-pass Butterworth filtered calcium fluorescence transients sampled at 1 kHz were used to identify a valid signal and reconstruct the AP train, respectively, as shown in Figs. 1 and 1. The number of zeros and poles in the Butterworth filter were set to be 3 and 10, respectively, at all times. The rising and decay duration triggered by a single AP firing usually last hundreds of milliseconds. Using the 10-Hz low-pass filter, the AP firing and corresponding rising durations are easily identified. The identified arising duration, the adjacent decay duration, and the following 500-ms duration consist of a part of the calcium fluorescence transients [Fig. 1], the time interval of which corresponds to the 100-Hz filtered signal. The corresponding signal is used to reconstruct the AP train, which is mostly based on the consideration that the rising duration trigged by two APs firing at a 25-ms interval can still be completely separated [Fig. 1] using a 100-Hz low-pass filter. Accordingly, MSF_TC can successfully identify APs firing at frequencies up to 40 Hz from the calcium fluorescence transients.

As indicated in the box in Fig. 1 and the middle panel of Fig. 1, the extracted rising duration p(t) t ∈ [t _i,t _{i+ 1+2k}] from the identified valid signal could meet the following criteria:

Eq. 1

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} p(t_{i + 1}) - p(t_i) > {\rm thr,} \end{equation}\end{document} $$p\left(t_{i+1}\right)-p\left(t_i\right)>\mathrm{thr},$$

Eq. 2

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{eqnarray} p(t_{i + 1 + 2m}) &<& \max [ {p({t_{i + 3 + 2m}}),p({t_{i + 5 + 2m}}),p({t_{i + 7 + 2m}})} ]\nonumber\\ m &= &0,\;1, \ldots, \;k - 1, \end{eqnarray}\end{document} $$\begin{array}{ccc}\hfill p\left(t_{i+1.2m}\right)& <& max\left[p\left(t_{i+3.2m}\right),p\left(t_{i+5.2m}\right),p\left(t_{i+7.2m}\right)\right]\hfill \\ \hfill m& =& 0,1,\dots ,k-1,\hfill \end{array}$$

Eq. 3

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} p(t_{i + 1 + 2k}) > \max \left[ {p\left({t_{i + 1 + 2k}} \right),p\left({t_{i + 5 + 2k}} \right),p\left({t_{i + 7 + 2k}} \right)} \right], \end{equation}\end{document} $$p\left(t_{i+1.2k}\right)>max\left[p\left(t_{i+1.2k}\right),p\left(t_{i+5.2k}\right),p\left(t_{i+7.2k}\right)\right],$$

where t _i is the time corresponding to the peak or valley value p(t _i), and thr denotes the threshold determined by clustering for the amplitudes of the noise and the calcium waves triggered by a single AP. A suitable setting is 0.6 to 0.7 times the amplitude of the template. The times [t ₁, t ₂,…, t _k], indicated by the arrows in the bottom panel of Fig. 1, correspond to the peaks of the differential of the extracted rising duration [Fig. 1, middle], and the peak values are larger than the standard deviation of the differential signal [horizontal line in Fig. 1, bottom panel]. These times are used to generate AP train candidates, as follows:

Eq. 4

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} \ell _s (t) = \left\{{\begin{array}{@{}l@{\quad}l} {1\;{\rm or}\;2} \hfill & {t = t_1} \hfill \\ {0\;{\rm or}\,1} \hfill & {t = t_i (i > 1)} \hfill \\ 0 \hfill & {{\rm otherwise}} \hfill \\ \end{array},} \right. \end{equation}\end{document} $${\ell }_s\left(t\right)=\left\{\begin{array}{cc}1\mathrm{or}2\hfill & t=t_1\hfill \\ 0\mathrm{or}1\hfill & t=t_i(i>1)\hfill \\ 0\hfill & \mathrm{otherwise}\hfill \end{array},\right.$$

where ℓ_s(t) is a candidate of AP train, s = 1, 2,…, 2^k. Here, ℓ_s(t _i) = 1 means that single AP fires at time t _i. Note that ℓ_s(t ₁) = 2 represents a special AP burst firing pattern that two AP fire at an about 10-ms time interval and their waves are not completely separated. This usually happens at the beginning of the burst.^{34, 35}

2.2.3.

Template convolution

The AP candidate train ℓ_s(t) was convolved with the template to get the corresponding simulating signal f _s(t). As shown in Figs. 1 and 1, the simulated signal [solid line in the Fig. 1] was obtained by convoluting the AP train candidate containing five APs [Fig. 1] with the template. The difference between the simulated signal f _s(t) and the extracted rising duration p(t) t ∈ [t _i, t _{i+ 1+2k}] is measured by the Euclidean distance shown in

Eq. 5

[TeX:] \documentclass[12pt]{minimal}\begin{document}\begin{equation} k_s = \mathop {\min}\limits_i \left\{{{\rm norm}\left[ {p_{[t,t + \Delta t]}, \;f_s (i,\;i + \Delta t)} \right]} \right\}. \end{equation}\end{document} $$k_s=\underset{i}{min}\left\{\mathrm{norm}\left[p_{[t,t+\Delta t]},f_s(i,i+\Delta t)\right]\right\}.$$

The minimum value k _s ^* indicates that the simulating signal f _s ^*(t) is the most similar to the rising duration. Thus, the AP train [TeX:] \documentclass[12pt]{minimal}\begin{document}$\ell _s ^ * (t)$ ${\ell }_s^{*}\left(t\right)$ corresponding to f _s ^*(t) is regarded as reconstructed AP train. As illustrated in Figs. 1 through 1, the AP train candidate in Fig. 1 is more suitable to represent the real AP train than that in Fig. 1.

3.

Results

The MSF_TC algorithm was applied to eight pieces of experimental data from four different cells. They are 300 s long and sampled at a 1-kHz frequency. The template of any trace was extracted from the individual data sets. The difference between the reconstructed and real AP trains reflecting the efficiency of the MSF_TC algorithm was weighed with false positive rates and event detection rates.

In Fig. 2, the data set 1 (top) contains 301 APs, and the visual comparison of real (middle) and reconstruction (bottom) AP train indicated that no difference between them exists. In this case, the false positive rate and the event detection rate were 6 to 11 and 94 to 95%, respectively, when the threshold was set to be 0.6 to 0.7 times the amplitude of the template, as indicted in Figs. 2 and 2 (labeled the number 1 on the x axis). As shown in Figs. 2 and 2, the reconstruction results from the other seven pieces of data were similar with that already considered.

Fig. 2

Capability of MSF_TC algorithm in AP train reconstruction: (a) calcium fluorescence transients (top), the corresponding real (Middle) and reconstructed (bottom) AP trains; (b) histogram of the difference between AP firing reconstruction and real timing obtained by analyzing data set in (a); (c) ISI distribution from the actual AP train [middle panel in (a)]; and (d) burst AP firing corresponding calcium fluorescence, real and reconstructed AP train. For an eight-data-set analysis, a 10% change of the threshold influenced by (e) the false positive rate and (f) the event detection. (g) plot of the standard deviation of the absolute error between real and reconstructed AP timings and (h) ISI distribution for eight real AP trains.

We also measured the time precision of the AP train reconstruction by comparing it with the real AP train. The absolute error of AP firing estimated timing and real timing distributed over the 0- to 10-ms interval was more than 90% [Fig. 2]. The standard deviation of the absolute errors was less than 7 ms for data set 1, and was also suitable for the other seven data sets [Fig. 2].

To show the ability of the MSF_TC algorithm to reconstruct the burst firing, a histogram of the inter-AP interval (ISI) less than 1000 ms was given by analyzing data set 1 [Fig. 2]. From the histogram, we can conclude that high-frequency and short-duration firing are common patterns in our data set. For further verification, ISI distribution from eight AP trains was also given and ISIs less than 50 ms held 30% [Fig. 2]. In addition, a part of calcium fluorescence transient triggered by a burst firing is shown in Fig. 2. During the AP burst firing, the rising duration generated by individual AP obviously changes, which heavily hinders the reconstruction of AP train using some characteristics of a calcium fluorescence transient, such as its rising amplitude and the amplitude variance. Considered from the reconstruction results and our data set features, the MSF_TC algorithm can effectively identify an individual AP from a burst.

To test the robustness of the MSF_TC algorithm against the threshold, a 10% change in the threshold was applied, namely, the threshold increasing from 0.6 to 0.7 times of the amplitude of the template. The false positive rate (FPR) and the event detection rate (EVR) were increased and decreased by an average of about 3%, respectively [Figs. 2 and 2].

We compared the reconstruction performance of the MSF_TC algorithm with that of two other typical methods, template matching^{22, 23} and the deconvolution method⁷ (Table 1). The calculations, including parameter selection, were carried out in accordance with previous work.^{7, 22, 23} For MSF_TC, the threshold value was set to 0.7 times the amplitude of the template. Low EDRs demonstrate that the template matching is unsuitable for reconstructing the neuronal burst firing. For the deconvolution method, the effective results are strongly dependent on the high SNR. For example, the SNRs of data sets 1 to 4 are about a range of from 2 to 2.5; in contrast, the data set 6 reaches 4. Therefore, compared with the deconvolution method, the MSF_TC algorithm could provide effective results from the calcium transients even with a relatively low SNR.

Table 1

MSF_TC algorithm compared with template matching and the deconvolution method.

4.

Discussion and Conclusions

As already noted, the neuronal bursts firings are extremely important for information encoding. We proposed a simple method to reconstruct an AP train from the calcium fluorescence transients that is suitable for even neuronal burst firing. By the multiscale filter application and valid calcium identification, the influence of long-duration noise on the reconstruction result is eliminated to a large extent. Furthermore, the simulation signals generated by the convolution of AP train candidates with the template can provide many patterns to match neuronal burst firing. These procedures guarantee the reliability and effectiveness of the MSF_TC algorithm in the analysis of the data sets.

Note that the method requires a proper template, so acquiring the template is an important procedure. Calcium fluorescence transients induced by a single AP is an all or nothing process.³⁶ In addition, the model for computing calcium dynamics demonstrates that the amplitude of the calcium fluorescence transient is almost constant.²¹ Therefore, the template can be successfully obtained by extracting and averaging the calcium traces generated by a single AP.

In the MSF_TC algorithm, multiscale filters are adopted. The scale selection is determined by the sum of the rise and decay duration formed by single AP firing and the noise type. There is a significant difference between the AP firing and no AP firing, corresponding to a 10-Hz low-pass-filtered calcium fluorescence transient. Thus, identifying valid calcium fluorescence is easy and the threshold is not discussed. In addition, properly increasing the length of the valid calcium fluorescence transient prevents missing the APs, whose firings lead to calcium wave part separations. Obviously, identifying an individual AP from a burst is impossible using the 10-Hz filtered signal. To resolve the rising durations triggered by two APs with 25-ms firing intervals, that is, the two rising durations are completely separate—100-Hz low-pass filters were properly selected, as verified by our results. Note that, in part of our data sets, the occurrence of two APs with an interval of 20 ms can cause the complete separation of the rising duration. As a result, the APs firing at a 50-Hz frequency from a burst were successfully identified.

Compared with other typical methods, i.e., template matching²² and the deconvolution method,⁷ the results illustrate that MSF_TC algorithm is superior to these two methods. The reason may be as follows: many of the burst firings corresponding to the calcium fluorescence transients lead the template matching to have a low efficacy and an outstanding performance of the deconvolution method is strongly dependent on a high SNR. For sparse AP firing, the MSF_TC algorithm is still inferior to template matching, because it chooses only the characteristics of the rising duration that identify an individual AP from a burst firing. Furthermore, like the other methods,^{9, 29, 30} the MSF_TC algorithm can not provide the AP firing rate when the rising durations are not separated. Properly use of decay duration may enhance the efficiency of identifying individual APs from a burst firing.

In summary, according to the characteristics of the data sets from our optical system, a reliable method was designed to reconstruct AP trains from calcium fluorescence transients. The results showed that the MSF_TC algorithm has a perfect performance in revealing an individual AP from calcium fluorescence transients corresponding to a burst firing, which should be beneficial for further deciphering neuronal circuits with a higher temporal resolution.