1. Introduction

With the development of fluorescent calcium indicators and fast imaging technology, calcium imaging has enabled us to optically detect neuronal activity within local neuronal populations [1–3], which is essential for understanding the neural information processing. Reconstruction methods could recover spike patterns from the recorded calcium traces, and play an important role in analyzing neuronal activity from optical imaging [4, 5]. Generally, detecting single spike firing from Ca²⁺ trace is an easy task because there are no overlapped Ca²⁺ waves. These non-overlapping Ca²⁺ waves can be easily identified with template matching. However, reconstruction of burst activity from Ca²⁺ trace is still a challenge. This challenge mainly originates from a large number of noise in Ca²⁺ imaging of neuronal population [6], high nonlinearities in the spike-calcium relationship [7–9], the often relatively low temporal resolution of the calcium signal when compared with the action potential timescale [10].

Significant progress has been achieved to address neuronal signal reconstruction [7, 8, 10–20]. The template matching method, which is considered a classical method for detecting weak signals, provides a consistent reconstruction of single-spike firing, but it is not suitable for burst firing reconstruction. Under the assumption that the recorded calcium fluorescence signal can be modeled by convolving the template calcium wave with the action potential (AP) train, the deconvolution (Deconv) method can reconstruct the burst firing from calcium fluorescence recording with high signal-to-noise ratio (SNR) [12]. Considering this assumption and some properties of the AP train like non negativity [8, 10] and sparsity, reconstructing the neuronal firings can be considered as equivalent to solving an optimization problem. This “conversion” of the problem to an optimization problem ensures that the reconstruction algorithm is suitable for analysis of signals with relatively low SNR [19]. Especially, a recently proposed method [16], which employs L1-minimization strategy, can automatically obtain AP train, the template and the baseline of signal from a given Ca²⁺ trace with relatively low SNR. Despite advances in neuronal burst activity reconstruction, most of the proposed methods incorporate only one or two features of spike trains in designing the optimization model for identifying spikes, which results in still no well-proven method suitable for complicated calcium traces with low SNR.

Here we propose a reconstruction method that applies Lq minimization approach and interval screening (LqIS). Lq minimization approach is a widely used technique for sparse signal reconstruction [21]. In the designed model for reconstruction, several properties of the AP train including non-negativity, sparsity, and firing rate are considered. The neuronal spike firings could be recovered from the complicated calcium fluorescence trace through Lq minimization. IS, introduced in our previous works [17, 22], is used to screen out the spike interval so that the reconstruction method can only analyze the calcium trace in these effective intervals. This procedure could largely reduce the dimension of the reconstruction model. Results show that the introduction of IS vastly increases the speed of reconstruction. Compared to some available methods, like L1, Deconv and non-negative deconvolution (NNDeconv), the LqIS is superior in reconstructing neuronal burst activity from Ca²⁺ traces with low SNR. Moreover, we demonstrate that spikes can be automatically reconstructed from the calcium imaging data of a neuronal population with the LqIS approach.

2. Methods

2.1 LqIS model for reconstruction of neuronal burst activity

In general, a calcium fluorescence trace, if no noise is presented, can be modeled as the convolution of a given template (corresponding to an AP) and a 0-1 sequence (AP train). Further, if the noise of the calcium fluorescence trace is considered to be close to Gaussian white noise, the calcium fluorescence trace can be expressed as the following equation [12].

here, ∆F/F (t) denotes the relative calcium fluorescence intensity change, P(t) the AP train, e(t) the noise, and Temp(t) the calcium waveform of an AP. By discretizing Eq. (1), we can obtain the following linear system

In this system, ∆F/F, P, and E are all time series sampled from (∆F/F)(t), P(t), and e(t) respectively. A_T is predetermined before reconstruction and is the matrix form of the time-shifted Temp(t) series in Eq. (1). Reconstruction of neuronal spiking activity, i.e., inference of the spiking timing from calcium fluorescence signals, is achieved by solving Eq. (2) under the condition of ∆F/F and A_T being known. However, Eq. (2) is ill-conditioned and cannot be directly solved. Considering the features of the spike timing, solving Eq. (2) can be considered as equivalent to solving the following optimization problem.

In this optimization problem, ||.||₂ is 2-norm and ||.||_q, q-norm, is defined as

P_i represents the i^th component of the AP train P, with n being the number of components in P; m represents the condition that a neuron can have up to one AP in the time interval of m milliseconds, and was set to 10 in this paper; λ is sparsity parameter, denoting the tradeoff between the fitting error and the sparsity of the AP train P. A large λ value will enforce the sparsity and increase the fitting error. And a smaller fitting error leads to a reduction in the sparsity to a certain extent. In the optimization problem represented by Eq. (3), there are two type of inequality constraints. One corresponds to the non-negative property of the AP train, and the other is used for constraining the highest spiking rate of a neuron, briefly named firing rate constraint. This optimization problem here is called the Lq minimization model. By solving this equation, we can determine a time series P that is close to the actual AP train.

In Eq. (3), we imposed Lq penalty on the coefficient, P, which is in a subspace determined by constraint conditions of Eq. (3). This increases the computational burden and requires a more complicated algorithm to solve it. To effectively solve Eq. (3), we designed an interval screening procedure to reduce the dimension of Eq. (3). Namely, we use screening procedure for estimating firing intervals whose union is denoted by M_. The optimization model (3) is rewritten as

where (A_T)_M is obtained by deleting the columns of A_T whose indexes are not in M, and P_(i∈M) represents a vector formed by P elements whose indexes are in M. The optimization model (5) is referred to as LqIS in this study. Through the IS procedure, the dimension of LqIS is reduced to the number of elements in M, which is far less than the dimension of Lq, i.e. the size of P. In comparison with solving Lq, solving LqIS can be more computationally economical. To reconstruct the AP train, we adopted previously proposed algorithms for solving the optimization problems (3) and (5) [21, 23]. Briefly, the optimization problem (3) or (5) was transferred into a series of convex optimization problem with the strategy presented by Simon Foucart, et al. [21], and each convex optimization problem was solved by the algorithm by Ziyou Gao, et al. [23]. Analysis of the effects of the parameter λ on the reconstructions shows when λ was set two in Eq. (5) the sparse reconstructions from the simulation data sets and experimental data sets could be obtained. In the following analysis, q in ||.||_q was set to 0.5, and sparsity parameter λ was set to 2.

We presented the procedure of IS for roughly estimating the time interval of neuronal spiking (Fig. 1). The procedure contains the following steps: (1) Smoothing the original calcium fluorescence trace using a low-pass filter and subsequently modifying the filtered trace by setting its negative value to zero; (2) Performing a differential operation on the modified signal and computing the peak area enclosed by the wave of the differential trace and the x-axis; (3) Selecting the waves whose peak area values are more than the given threshold value, that is, 0.15 times the amplitude of the template; and (4) Generating the interval candidates by thresholding AP spiking score. For each selected wave in step 3, we regarded its corresponding time interval in which the signal value is higher than the threshold as the estimated interval of AP spiking, as indicated by the arrows in Fig. 1. Then AP spiking score of each selected wave, which is an index to depict the possibility of AP spiking in the estimated interval, is obtained by solving the following optimization problem:

 

Fig. 1 Steps of LqIS procedure for neuronal burst activity reconstruction from calcium fluorescence trace. (1) Ca²⁺ is modeled as the sum of the convolution signal and Gaussian noise. (2) The smoothed trace is obtained with a 100-Hz Butterworth filter. (3) The peak areas in the differential modified trace are calculated and then used to screen the AP interval. (4) For each candidate spiking interval obtained from the screening procedure, the spiking score (inverted triangle) is calculated by solving the non-negative deconvolution model. (5) By thresholding the spiking scores, the AP spiking intervals used for reconstructions are selected. (6) From the determined AP spiking intervals, the AP firing timings are calculated by solving the Eq. (6).

Download Full Size | PDF

here, (A_T)_m and P_(i∈m) represent the corresponding parameters in optimization problem (5), while m represents the sampling time points in the estimated interval. In our study, the optimization problem (6) was solved by using the gradient descent algorithm [24]. The AP spiking score of interval m is defined as the sum of all the elements of P_m, as indicated by the red inverted triangles in Fig. 1. The time interval over which the AP spiking score is more than 0.35 is considered as the screened time interval that is used in the optimization problem (5), and here, the threshold value is experimentally determined as the value that ensures that almost all the spikes fall in the range of the screened time interval. Figure 1 shows an example of the screened time intervals. Note that, in designing IS procedure, we considered two features of Ca²⁺ traces: one is that AP spiking appears the rising duration of Ca²⁺ trace [12], the other is some notable differences between Ca²⁺ waveform of single AP and noise.

To illustrate that LqIS has a sparser reconstruction than L1 minimization method [16], we simulated traces with a fixed template and a signal baseline for this purpose. In this case, we described L1 minimization method as follows

The symbols in the above optimization problem have the same meanings as those in Eq. (4). We employed a sparse reconstruction algorithm [24] to solve Eq. (7), and the sparsity parameter λ was experimentally determined and set to 2.

2.2 Selection of two parameters and template estimation

In LqIS method, two key parameters are toned to be predetermined. One parameter is spiking score threshold, and the other is sparsity parameter in Optimization problem (3) and (5). We generated the simulation signals corresponding to a single AP spiking, calculated the spiking score of the interval where AP appeared, and depicted the frequency distribution of the calculated spiking scores (Fig. 2(a)). We noted that the distributed range of spiking scores increase as the SNR decreases. For low SNR (SNR = 1.5), the minimum value among spiking scores was about 0.65, far more than our set value, 0.35. In the case of neuronal burst activity, the spiking scores in the burst interval can be contributed by more than on AP firings and will be even larger (larger than 0.80, as simulated with two merged spikes shown in Fig. 2(b)). In the procedure of IS, the threshold of spiking score was used to screen out the intervals for AP train reconstruction. The above results indicated that our set value, 0.35, can screen out all the time intervals containing AP spiking. At the same time, setting a small value means that the spiking score setting would generated some screened time intervals without AP spiking inside. This case is permissible, as shown in 5th panel of Fig. 1. As for the sparsity parameter, results in Fig. 2(c) illustrated the strong robustness of sparsity parameter on the reconstructions. We noted that the reconstructions only slightly changed when sparsity parameter lies in the range from 1.25 to 8, when SNR is set to 2 and 5.

 

Fig. 2 Selection of two parameters in LqIS. (a) shows the distribution of spiking scores of single spike traces. The simulation traces were classified into four groups based on their SNRs. Each group contains 200 traces with the same SNR, each trace contains one single spike. (b) The distribution of spike scores of the traces with two merged spikes. (c) The reconstruction value, represented by the grayscale lines, varies with the sparsity parameter λ. Top, a simulated trace generated by the convolution of the AP train containing three APs with a natural exponential function. Middle, reconstruction values from the top simulation trace when noise was added (SNR = 2). Each row represents values from one trace with one given sparsity parameter λ, each column shows reconstruction values of the same signal with different sparsity parameters. Bottom, reconstruction values from the top simulation trace with high SNR (SNR = 5).

Download Full Size | PDF

Template estimation is also a key factor to AP train reconstructions. LqIS reconstructs AP train using a fixed template. In our simulation analysis, the template was the one used for generating the data set. In the experimental data set analysis, the template was the average of 10 isolated Ca²⁺ waveforms, all of which correspond to a single AP spiking. We extracted Ca²⁺ waveforms with manual operation.

2.3 Generation of simulation data sets

Simulation traces are generated based on the assumption that the calcium fluorescence signal can be acquired by the convolution of the calcium fluorescence template with an AP train [12, 25]. Here, the template is a natural exponential function with a decay time of 100 ms, and the AP train is a 0-1 series of 1000 ms long. In the AP train, the timings corresponding to a value of 1 are considered as the spike timings. A set of simulation calcium traces were generated as the sum of Gaussian noise and the convolution of the template and the AP train. The SNR of the calcium trace is the ratio of the template amplitude to the standard deviation of the Gaussian noise.

2.4 Pretreatment of experimental data sets

Before the reconstruction, some pretreatments were performed on done to the experimental data sets. Firstly, the calcium traces were smoothed with a low-pass Butterworth filter (100 Hz, with the number of zeros and poles set to 3 and 10, respectively). Secondly, from the 100 Hz filtered signal, the waves corresponding to neuronal firing were extracted. The template was calculated by averaging these extracted waves, since it is required when using LqIS and non-negative deconvolution (NNDeconv) for the reconstruction. Thirdly, considering that solving the optimization problems in LqIS and NNDeconv are time-consuming for calcium traces with large data sizes, we estimated the time interval of neuronal firing from the given signal, and we used LqIS and NNDeconv to analyze only the 100-Hz filtered signals of the estimated interval. The estimating procedure has been proved to be effective in our previous study [17], and here, it is briefly described. An extremely low-pass Butterworth filter (10 Hz, with the number of zeros and poles set to 3 and 10, respectively) was used for smoothing the original calcium fluorescence signals. From the 10-Hz filtered signal, the rising duration [t₀, t₁] was extracted via the threshold method [17, 22], with t₀ and t1 denoting the start and end points, respectively, of the rising duration; here, the interval [t₀, t₁ + 200 ms] was regarded as the estimated interval of neuronal firing.

2.5 Quantification of reconstructions

By comparing the AP train reconstructed from the calcium fluorescence trace with the AP train extracted from the neuronal electrical signal by using the threshold method, we evaluated the reconstruction methods using recall rate and precision rate. Here, the recall rate is defined as the ratio of true positive APs from the detected APs to the true APs observed in the neuronal electrical signal. The precision rate is the ratio of the true positive APs to the recognized APs. A recognized AP is regarded as a true positive AP when the time difference between the recognized and true positions is less than 10 ms for a calcium fluorescence trace with more than 100-Hz sampling frequency, and 15 ms for a calcium fluorescence trace with a 100-Hz sampling frequency.

3. Results

3.1 Validating the procedure of IS

To validate IS in the reconstruction of neuronal burst activity, calcium fluorescence signals of 1 s were generated by convolving the calcium fluorescence template with the AP trains. To address every possible scenario, in this simulation, the AP train was set to contain 6 action potentials with firing intervals of 20, 30, 40, 50, and 260 ms respectively, and the SNR of the simulated calcium fluorescence signal was set to 3.5. We generated four blocks of calcium fluorescence signals corresponding to sampling frequencies of 1000, 500, 200, and 100 Hz. Each block contained 20 calcium fluorescence signals. Figure 3(a) shows the calcium fluorescence signal sampled at 1000 Hz and the APs reconstructed with IS (LqIS, Eq. (5) and without IS (Lq, Eq. (3)). From these two reconstructed AP trains, we observe that the LqIS method does not degrade the reconstruction accuracy. This conclusion was also further verified by using LqIS and Lq to analyze the four blocks of calcium fluorescence signals. The quantitative results are shown in Figs. 3(b), 3(c), and 3(d). In Fig. 3(b), for each block of calcium fluorescence trace, the average recall rates of the reconstructions derived from LqIS (red dots) and Lq (black dots) are depicted. Kolmogorov-Smirnov Test results show there is no statistical difference between reconstruction recall, precision and timing error between with and without IS procedure. The recovery speed of LqIS is 20- to 50-fold faster than that of Lq, as can be observed from Fig. 3(e). So we can conclude that LqIS can quickly reconstruct neuronal burst activity from calcium fluorescence traces without compromising on reconstruction accuracy. Here, the computational time of LqIS is the sum of running time of all steps, including the IS procedure (Steps 3-5 in Fig. 1). Considering that simulated data sets were mainly made of single spike firing (Fig. 3(a)), and the recall rate and precision of reconstruction on these data sets from Lq and LqIS showed little difference (Fig. 3(b)-3(d)), it is concluded that IS procedure would not miss single spike firing.

 

Fig. 3 Introduction of IS enhances the reconstruction speed without compromising the reconstruction accuracy. (a) The performance comparison of Lq minimization with and without IS shows that the introduction of IS has no appreciable effect on the reconstruction. (b)-(c) Means and variances of recall and precision rate derived from LqIS and Lq, respectively. (d) The mean and variance of timing error by LqIS and Lq. (e) The total computational time of 20 simulated calcium traces reconstruction with LqIS and Lq.

Download Full Size | PDF

3.2 Strong sparsity of reconstructions with LqIS

To verify that LqIS has the ability to sparsely reconstruct from Ca²⁺ signals, we compared the reconstruction results from LqIS and L1 minimization on the simulated traces sampled at 500 Hz. Figure 4 shows two reconstruction examples (Fig. 4(a)) and the reconstructed signals on 200 sets of Ca²⁺ signals (Fig. 4(b)). In the reconstructions, we used a large value of λ in Eq. (5) and Eq. (7) for obtaining the sparse reconstructed signal. Results shows that for traces with low SNR (SNR = 2), L1 minimization method just obtains weak sparsity of reconstructions (Fig. 4(b)), gaps still exist between results from L1 and the real AP train. Furthermore, with L1 minimization, reconstruction value of resting fluorescence period has similar signal intensity with that of AP firing, which would be regarded as real spikes. LqIS sparsely reconstructed the AP train even from low SNR traces (Fig. 4(b) top), and the reconstruction agreed with the real AP train. In high SNR cases, Lq provided reconstructions nearly equal to the actual AP train.

 

Fig. 4 Comparisons of the reconstructions derived from LqSIS and L1 minimization. (a) Two simulated traces and the normalize spike trains from LqIS and L1 minimization methods. (b) Reconstructions derived by LqIS (top) and L1 minimization (bottom) on 200 trials. Each row represents a piece of reconstructed signals and signal intensities are displayed in grey level.

Download Full Size | PDF

In addition to Lq norm penalty, firing rate constraint on the AP train can also contribute to the reconstruction sparsity. We compared the reconstructed signals from LqIS to that from Eq. (3) without firing rates constraint, on the test signals with SNR = 2 and sampled at 500 Hz (Fig. 5). Without firing rate constraint, the reconstructed signal intensities are sometimes greater than one for better approximating the original signal, i.e., to achieve a smaller least-square error in Eq. (3), as shown in Fig. 5(a). This is a common phenomenon (Fig. 5(b)). With firing rate constraint, the sparsity of the reconstructed signals was obviously enhanced. The reasons may be originated from the fact that firing rate constraint narrowed the range of the solution of Eq. (3), which may make Lq norm penalty more effective, in the meantime this constraint was consistent with the feature of AP train.

 

Fig. 5 Necessity of firing rate constraint in LqIS. (a) A simulation trace and their corresponding recovery when with and without firing rate constraint. (b) Use the data sets, containing 200 sets of Ca²⁺ simulation signals, to verify that firing rate constraint can enhance the sparsity of the reconstructed signal. For each set of Ca²⁺ simulation signal, two reconstructions derived by LqIS with and without firing rate constraint (row) display. The intensities of the reconstructed signal are displayed in grey level.

Download Full Size | PDF

3.3 Performance of LqIS in reconstruction of simulation data

We quantified the reconstruction results obtained with LqIS from simulated calcium traces with SNRs ranging from 1 to 3. For a fixed SNR and firing rate, a corresponding block of signals was generated that contained 50 sets of calcium fluorescence signals with an interval of 1 s and sampling frequency of 1000 Hz. The AP trains and the corresponding calcium fluorescence signals are shown in Fig. 6(a). From Fig. 6(a), we observe that when the SNR is 2, manual detection of the AP firing from the calcium fluorescence signals is difficult. However, from Figs. 6(b)-6(d), we note that when the SNR is 1.5, LqIS can generate perfect reconstructions (recall and precision of about 0.95) for different firing rates; when the SNR is below 1, LqIS still provides effective reconstructions. These results indicate that LqIS is robust to noise when reconstructing AP burst firing along with providing millisecond precision. We also point out that in an AP firing range between 15 Hz and 35 Hz, the reconstructions derived from LqIS are very similar (Figs. 6(b)-6(d)), probably because the high sampling frequency (1000 Hz) provides sufficient information for the reconstruction of the high AP firing rate.

 

Fig. 6 Neuronal activity reconstruction with LqIS shows that LqIS is robust to noise and large spiking rates. (a) Action potential (AP) trains with various spiking rates and their simulated calcium fluorescence signals. (b-d) The means and variances of recall, precision, and timing error as obtained using LqIS from calcium traces with varying signal-to-noise ratios (SNRs) and AP spiking rates.

Download Full Size | PDF

We also quantified the effect of the data sampling frequency on the reconstruction by analyzing the simulated traces via LqIS. The simulated calcium trace is modeled as the convolution of the template with AP trains with different intervals, as shown in Fig. 7(a). The convolution signals with noise were sampled at frequencies of 1000 Hz, 500 Hz, 200 Hz, and 100 Hz, and the four corresponding sets of signals are shown in Fig. 7(b). This approach allowed us to generate blocks of calcium fluorescence signals each containing 50 sets, with each set having the same sampling frequency and SNR. For each block of signals, we calculated the recall, precision, and timing error of the reconstructions derived from LqIS. The results shown in Figs. 7(c)-7(d) indicate that the reconstruction accuracies of neuronal burst activity decrease as the sampling frequency decreases, mainly because low-frequency sampling deteriorates the signal information of neuronal burst firing. Despite this deterioration, LqIS still provides effective reconstructions in the high-SNR cases (SNR > 2.5) when the sampling frequency is below 200 Hz, as shown in Figs. 7(c) and 7(d).

 

Fig. 7 Effect of the sampling rates of calcium fluorescence signals on the reconstruction accuracy. (a) Action potential (AP) train and its corresponding simulated calcium fluorescence signal without noise. (b) Simulated calcium trace sampling at 1000 Hz, 500 Hz, 200 Hz, and 100 Hz. (c)-(e) Means and variances of recall rate, precision, and timing error of the reconstructions from calcium fluorescence signals with different sampling frequencies at different SNRs.

Download Full Size | PDF

To illustrate the merits of LqIS, we analyzed simulated calcium traces with different SNRs that were obtained using the LqIS, L1, NNDeconv, and Deconv approaches. In the simulation, the AP trains have the same settings as in Fig. 3, and for a fixed SNR, 50 sets of calcium fluorescence signals with time intervals of 1s and sampling frequency of 1000 Hz were used for our quantification analysis. Figure 8(a) illustrates two sets of calcium fluorescence signals with SNR = 1 and SNR = 3 along with the corresponding reconstructed signals derived from the abovementioned methods. From Fig. 8(a), we observe that for calcium fluorescence signals with high SNR, all four algorithms generate perfect reconstructions. However, in the low-SNR case (SNR = 1), Deconv, NNDeconv and L1 generate reconstructed signals from which identifying the AP firing timing is difficult; in contrast, LqIS provides a reconstructed signal that faithfully reflects the AP train. This ability to effectively reconstruct the AP train from low-SNR calcium traces forms one of the key advantages of LqIS. Further statistical tests (Kolmogorov-Smirnov Test) show that, LqIS has a significant advantage (p < 0.005) in recall rate compared to other three methods at low SNR (SNR = 1.0, 1.5 and 2.0) in Fig. 8(b).

 

Fig. 8 Performance of the four reconstruction methods with simulated signals with signal-to-noise ratios (SNRs) of 1 and 3. (a) Two sets of simulated calcium fluorescence signals with SNRs of 1 and 3, and their normalized spike trains from LqIS, L1, NNDeconv and deconvolution (Deconv). (b)-(d) Means and variances of recall, precision, and timing error of the reconstructions from 50 data for quantifying the abovementioned methods.

Download Full Size | PDF

3.4 Performance of LqIS in reconstruction of experimental data sets

We also applied LqIS to experimental data sets to demonstrate the merits of the proposed algorithm. In this study, acute cortical brain slices from postnatal (P14−20) Wistar rats were used for the experiments. To induce epileptiform discharging, a potassium channel blocker 4-aminopyridine and high potassium solution were applied to the bath at final concentrations of 100 μM and 10 mM, respectively [26]. To compare the reconstructed with the real spike train, neuronal Ca²⁺ transients and spike firing were simultaneously recorded. Spike firing, sometimes referred to as membrane voltages, was recorded with a patch-clamp amplifier (EPC-9, HEKA). The recording pipettes were routinely filled with electrode filling solution containing 100 µM Fluo-5F (calcium dye). The local neuron population was loaded with the calcium fluorescence dye OGB-1/488 AM with the bolus-loading approach [2]. In this study, a custom-constructed fast two-photon fluorescence microscope system [27] was used, in which an AOD (Acousto-Optical Deflector) [28, 29] enables us to fast random scanning and recording Ca²⁺ fluorescence transients from the interested regions. The imaging process was controlled by custom-written software written in LabVIEW. Relative fluorescence intensity fluctuation (ΔF/F) was used to indicate neuronal firing.

Five sets of 300-s calcium fluorescence traces sampled at 1000 Hz were processed with this method, and the results are shown in Fig. 9. Figure 9(a) shows one period of a calcium trace and its corresponding electrical signals. Figure 9(b) illustrates three typical calcium transients corresponding to single-AP firing, two merged APs firing, and AP burst firing in order to verify the efficacy of reconstruction. From these reconstructions, we observe that the Deconv method poorly reconstructs AP firing at millisecond-level precision. Although the NNDeconv method performs better than the Deconv method, the extremely dense AP firing and calcium signals severely contaminated by noise are still difficult to be recovered. L1 method has a similar reconstruction as NNDeconv. These disadvantages of NNDeconv can be overcome to some extent by LqIS due to its constraints of the strong sparsity of the AP train and AP firing rates, and the LqIS results indicate that LqIS can provide a better reconstruction than NNDeconv. This conclusion is further validated by quantifying the reconstructions with the four methods, as shown in Figs. 9(c)-9(e). The reconstructions are from five data sets that contain 128, 110, 53, 87, and 300 APs in that order, as detected from recorded electrical signals. From Fig. 9(c), we observe that the recall rate for LqIS ranges from 0.91 to 1.0, which is higher than the recall rate for the NNDeconv method (0.7 to 0.96). From Figs. 9(d)-9(e), we note that LqIS, NNDeconv and L1 exhibit high precision rates and low timing errors. From these quantitative comparisons, we conclude that LqIS produces better reconstructions than the other three methods. Through the comparison among these methods, the advantage of LqIS is more obvious in the reconstruction of experimental data sets than the simulated data sets (Fig. 8). We considered that this case might be caused by the high frequency or nonlinear AP firings: the experimental data sets include high frequency or nonlinear AP firings, which results in the Ca²⁺ waves not being well described with a convolution model, in this case, LqIS considers more features of AP train than L1 and NNDeconv, and thus has a stronger ability to reconstruct the AP trains.

 

Fig. 9 Performance of various reconstruction methods with experimental signals. (a) Calcium fluorescence signal and its corresponding electrical signal. (b) Signals processed using LqIS, L1, NNDeconv, and Deconv. From top to bottom: calcium trace, corresponding cellular voltage, and their normalized spike trains from LqSIS, L1, NNDeconv, and Deconv methods, respectively. From left to right: three calcium traces corresponding to a single spike, two merged APs, and a string of APs. (c-e) Recall, precision, and timing error of reconstruction results from five experimental data sets as obtained with the above methods, respectively.

Download Full Size | PDF

Finally, we examine the LqIS performance in the reconstruction of complex neuronal population activity from calcium traces (Fig. 10). Figure 10(a) shows the fluorescence image of the population labeled with the calcium dye OGB-1/488 AM, and the simultaneously recorded traces of 50 neurons are displayed in Fig. 10(b). Each line of these 50 traces has a time interval length of 140 seconds with a sampling frequency of 1000 Hz; thus, these 50 traces totally contain 7 million data points. Figure 10(c) shows the reconstruction of the AP firing trains of the neuronal population obtained using LqIS without any human interference. From this figure, we observe that the neuronal population burst activities exhibit a synchronization similar to that of the calcium fluorescence traces in Fig. 10(b). For a better visualization, a magnification of the green box in Fig. 10(c) is shown in Fig. 10(d). For the reconstruction shown in Fig. 10, the time required by LqIS is about 4 h (with an Intel(R) Core(TM) i7-3632QM CPU@ 2.20 GHz). It is to be noted that this time estimate is based on MATLAB programing; on the other hand, use of the C\C++ code for LqIS can significantly decrease the calculation time. From these results, we can conclude that LqIS is suitable for neuronal population activity reconstruction.

 

Fig. 10 Reconstruction of neuronal population with LqIS. (a) Fluorescence image of population labeled with calcium dye OGB-488/AM. (b) Calcium fluorescence traces of 50 neurons are displayed to show oscillation in the population. (c) The reconstructed action potential (AP) firing from the population calcium traces. (d) Magnification of the trains shown in the green box in (c). The scale bar in (a) is 20 μm.

Download Full Size | PDF

4. Discussion and conclusion

In this study, we proposed the LqIS method for reconstructing the neuronal burst activity from calcium fluorescence traces of a neuronal population. Our reconstruction results show that the LqIS method can suitably reconstruct neuronal activity from calcium fluorescence traces with low SNRs and large data size. In comparison with the L1, NNDeconv and Deconv methods, the LqIS improves the accuracy of reconstruction. These results suggest that LqIS can be utilized for inferring the firing pattern of a large-scale population of neurons from their recorded calcium fluorescence signals.

Given that calcium fluorescence can be simply acquired by convolving a fixed template with the corresponding action potential train, the neuronal burst activity can be derived from the calcium fluorescence traces, essentially by solving a linear system of linear equations. However, it is difficult to directly solve the linear system due to its ill-conditioned property. When combined with certain features of the AP train, such as non-negativity, sparsity, and the firing interval, the linear system of equations can be converted into an optimization problem. It has been proved that this conversion boosts the reconstruction accuracy [19]. Unlike previous methods in which only one feature of the AP train is considered, our proposed method employs the features of non-negativity, sparsity, and the firing interval, which result in a more sophisticated optimization problem. LqIS employed more features of the AP train, so this method can generate reconstructions that are closer to the AP train.

Different from commonly L1 minimization method, our LqIS carried out the Lq penalty on the coefficients with the solution space being effectively narrowed by more constraints. These constraints were well consistent with the features of AP train, and thus the reduction of solution space was reasonable. This strategy assured that LqIS could provide the sparser reconstructions than L1 minimization method especially for the Ca²⁺ signals with low SNR (Fig. 4). This strategy obeys a simple rule that reducing solution space based on the features of the solution can really have more chance for searching the real solution. Note that this strategy has its disadvantages that the AP reconstruction becomes solving a sophisticated optimization problem and many fast algorithms for L1 minimization problem can’t be directly used.

The introduction of IS procedure increases the reconstruction speed of neuronal burst activity and is based on two considerations. On one hand, as described above, to solve a sophisticated optimization problem usually is very time consuming. The introduction of IS could reduce the dimensions of the optimization problem and thus boost the speed of the solution, as indicated in Fig. 3(e). On the other hand, the advances in calcium fluorescence labeling and imaging techniques enable us to record the calcium signals of a large neuronal population with millisecond resolution, which would produce a certain complicated data set with large data size. In this case, it is impractical to solve the introduced optimization model without IS. Thus, it becomes necessary to apply IS to reduce the dimensions of the optimization model.

The difference between calcium transients and shot-noise enables us to roughly estimate the firing durations with the use of a low-pass Butterworth filter, which is valuable for employing the sophisticated optimization problem to reconstruct the neuronal spiking train. The proposed IS was designed for this identification, and it was proved to be very effective. This IS strategy can be also applied in other reconstruction methods to boost the reconstruction speed. With increase in the size of the imaging data, the demand for large-sized data processing has also increased; in this context, IS can find wide applicability.

Estimating the algorithm error is also an important consideration in neuronal population burst reconstruction. In our previous work, we quantified the reconstruction errors derived from the local differential method by generating a simulated data set whose features were close to actual traces. However, the local differential method is inferior to LqIS in analyzing calcium fluorescence traces with low signal-to-noise ratios and long durations of burst activity (data not shown). Thus, in the context of our previous work [22], LqIS is suitable for the reconstruction of neuronal population burst activity.

We also introduced IS procedure into L1 and NNDeconv methods, named L1IS and NNDIS respectively, to further demonstrate the function of IS in for the reconstruction of burst activity. We used the same spike train and template as shown in Fig. 3(a) to generate the simulated Ca²⁺ traces with sampling frequency of 1000 Hz, each data set includes 20 Ca²⁺ traces. An example is shown in Fig. 11(a). From the Ca²⁺ trace in Fig. 11(a), the outputs of L1 and L1IS are almost the same (Fig. 11(a)). The similar results are also obtained for NND and NNDIS. We quantified the difference between the outputs of L1 (NND) with and without IS procedure. We found that this difference is negligible at different SNRs (Figs. 11(b)-11(c)). Comparison of time consumption showed that L1IS (NNDIS) enhanced the reconstruction speed more than 3 times in comparison with L1 minimization (NND), as shown in Fig. 11(d). So we could draw a conclusion that IS procedure do boost the reconstruction speed without any sacrifice of reconstruction accuracy. When comparing the recovery speed among these methods, we found that L1IS and NNDIS are about 6 times faster than LqIS (Fig. 3(e) and Fig. 11(d)).

 

Fig. 11 Introduction of IS enhances the reconstruction speed of L1 and NNDeconv methods. (a) The simulated Ca²⁺ trace and its corresponding AP train, and recovering spike train through L1, L1IS, NND and NNDIS methods. (b) The difference between the outputs of L1 (NND) with and without IS procedure is expressed as root mean square deviation, and the deviation refers to the Euclidean distance between two outputs. Each data set includes 20 Ca²⁺ traces. (c) Difference described in (b) is measured with Pearson correlation coefficient. (d) The total computational time of L1 and NND when with and without IS procedure. IS enhances the reconstruction speed over 3 times in L1 minimization and NNDeconv methods.

Download Full Size | PDF

In conclusion, the LqIS method yields both increased accuracy and speed of AP train reconstruction. The introduction of the Lq minimization model ensures immunity to noise, leading to enhanced recall accuracy. In addition, the introduction of IS reduces the dimension of the Lq minimization model, which significantly improves the reconstruction speed. Consequently, the method is applicable for calcium fluorescence signals of large sizes. These two advantages of the LqIS method ensure effective reconstruction of the electrical activity of a large-scale neuronal population, which might be helpful to the study of neuronal circuits.