2.1

Spiking sorting by Bayesian statistical model

In this section, the results of algorithm for traditional waveform classification are taken as prior knowledge. and Bayesian statistical model is used to optimize the classification results, identify the variable waveforms that are missed or falsely detected, and accurately obtain the spiking trains of single neurons. According to the results of the algorithms for traditional waveform classification, the prior probability function of the spiking trains d_j is:

where d_j(t) is a binary variable. When neuron j evokes spikes at time t, d_j(t) = 1, otherwise, d_j(t) = 0.

Assuming that the membrane voltage recorded by the electrode at the time t is a linear superposition of all the spiking waveforms in the neuron cluster by different forms at time t. The generative model of the membrane voltage at time t is established as follow^[7]:

where, the column vector represents the spikes recorded by all electrodes at time t, n_c is the number of muti-electrodes, and represents the spikes recorded by i th electrode at time t. All spiking waveforms are discretized into n_τ parts on average, and the matrix denotes the τ th segment waveform of neuron j recorded by all electrodes. is the background white noise with zero mean.

According to (1), the generative model for spikes throughout the entire observation period is obtained:

where, W * D represents the convolution operation in (2), and η is a multivariable Gaussian white noise. Therefore, V − W * D follows a multivariate Gaussian distribution. Consequently, under the conditions of spike times and spiking waveforms, the conditional probability function of the spikes generated is:

where Λ is the covariance matrix of the multidimensional Gaussian variable η, which not only represents the covariance of noise over time but also the covariance of noise caused by the spatial positions of the electrodes in the collection process.

Based on the prior function, the conditional probability function, and Bayesian theory, the posterior probability function of spikes evoked by the spiking trains and spiking waveforms is:

The spiking trains and the spiking waveforms are estimated by making the posterior probability obtain its maximum value[16-18].

2.2

Information coding by probabilistic statistical model

The spiking activity of neurons is affected by many factors at the same time, and the study of neural activity usually involves three different covariates. First, neural spiking activity is often associated with external biological and behavioral covariates, such as sensory stimuli and behaviors or specific motor outputs[13]. Second, the neuron’s current spiking activity is also related to its own spike-history[19]. Finally, the advancement of modern technology allows simultaneous record spiking activity of multi-neuron, and show that the spiking activity of a particular neuron is related to the spike-history of other neurons[15].

In this paper, a probabilistic statistical model is used to study the encoding and decoding process of acupuncture neurons. Firstly, the encoding model of acupuncture neuronal response was established using the framework of linear-nonlinear-Poisson process. In linear filtering, the external covariates, spike-history of self and the spike-history of other neurons are analyzed simultaneously. The three components undergo linear filtering and are then subjected to a nonlinear exponential filtering to generate the intensity function of the Poisson process[20,21]:

where k_i is the coupling parameter of external stimulation to neuron i, h_ii is the coupling parameter of own spike-history of the neuron i, and h_ij is the coupling parameter of other neurons. Let θ_i = (k_i, h_ii, h_ij), is the parameter set of the encoding model.

Assuming that the number of spiking events occurring in the time window dt is n_t, then there is n_t ~ Poiss[λdt]. Thus, the logarithmic probability function generated by neuronal spiking trains is as follows:

where c is a constant that does not depend on the parameter θ.

Then, the MLE was used to fit the spiking trains evoked by acupuncture to the model and estimate the model parameters. The logarithmic probability function log P(D | θ) is the logarithmic likelihood function of the model, and we can obtain the MLE of the model parameters by setting it to its maximum value[22,23].

3.1

Classification of acupuncture electrical signals

Using the spiking activity of neurons in the dorsal horn of the spinal cord from an acupuncture experiment on the Zusanli (ST36) acupoint in rats as the raw data. We perform both traditional classification algorithms and Bayesian statistical model based spiking classification, and compare their classification effects. In the acupuncture experiment, healthy adult rats were selected as experimental subjects, the acupuncture frequencies were 30 times/min, 60 times/min, 120 times/min and 180 times/min respectively.

Taking the stimulation frequency of 60 times/minute as an example, we selected the classification results of the two algorithms in 0.3 seconds to compare, as shown in Figure 1. Fig. 1 (a) shows the spiking trains of two neurons detected by the traditional classification algorithm. Fig. 1 (b) shows the spiking trains of two neurons detected in the Bayesian statistical model. Fig. 1 (c) shows the raw data collected by the electrodes. In the figure, the blue dots are the spikes extracted by both algorithms, the red dots are the spikes extracted by the Bayesian statistical model, and the green dots are the spikes omitted by the Bayesian statistical model.

Figure 1.

Comparison of the results between the two classification algorithms.

From the Fig. 1, it can be seen that the Bayesian statistical model has effectively solved the problem of superposition of multiple spiking waveforms. In particular, in the red dotted box, the spikes recorded by the electrode are concentrated. It shows that the responding activity of neurons is intense, and the Bayesian statistical model can detect all the spikes. In addition, when dealing with noise, such as the spiking waveforms corresponding to the spike in the green dotted line box is within the background noise range, such spikes are eliminated. It can be seen that the Bayesian statistical model reduces the missing rate of traditional classification algorithms and greatly improves the accuracy of classification results.

Under the conditions of acupuncture at four frequencies, the number of spiking events of two neurons estimated by the two classification algorithms is shown in Table 1. It can be concluded that the number of spikes detected by these two neurons in the Bayesian statistical model is significantly increased, and some even twice the number of spiking events in the traditional classification algorithm.

Table 1.

The number of spiking events in two classification algorithms.

In addition, we found that after the stimulation was increased to 120 times/min, although the frequency of the stimulation continued to increase, the production of spikes did not increase anymore. This is due to the limited number and availability of ion channels on neuron cell membranes, as well as the saturation of neurotransmitter release. When the frequency of stimulation is too high, the release of ion channels and neurotransmitters reaches saturation, resulting in no more spikes being generated. This phenomenon helps to maintain the stability and accuracy of the nervous system, ensuring the reliability of nerve signals during transmission. In medical research and practice, the stimulation frequency of neurons can be adjusted according to this principle to achieve the purpose of treatment and research. In clinical treatment, the commonly used frequency of stimulation is 30 times/minute to 120 times/minute. Therefore, in the follow-up study, we mainly analyzed the neuronal spiking activity under the three stimulus frequencies of 30 times/minute, 60 times/minute and 120 times/minute.

3.2

Encoding and decoding of acupuncture spiking activity

We verify the accuracy of the MLE through a set of simulation data. First, the intensity function was used to generate simulated spiking trains, and the model parameters were set to the real values in Table 2. Then, the MLE is used to fit the model and estimate the parameters of the model. The estimated results are shown in Table 2. By comparing the real value and the estimated value in Table 2, the relative error distribution of the estimated model parameters ranges from 0.62% to 4.47%. It is concluded that the MLE can accurately fit the model and estimate the model parameters.

Table 2.

Real and estimated values of model parameters.

In order to further study the difference of neuronal responding under the condition of acupuncture manipulation with different frequencies, a probabilistic statistical model was used to encode and decode the neuronal responding activity evoked by acupuncture of different frequencies, and the coupling parameters were estimated by MLE. The estimated results are shown in Table 3. We found that the coupling parameters of stimulus k were significantly lower than the coupling parameters of spike-history h under all three frequencies of acupuncture. With the increase of frequency, the more historical spikes of neurons, the smaller the coupling parameters of stimulus. It is concluded that the spike-history is the main factor to evoke the neuronal response in the experiment of neuronal spiking activity of acupuncture.

Table 3.

The coupling parameters of different acupuncture frequencies.