Method for prediction of cortical spiking trains

Abstract

The present invention is related to a method for prediction of cortical spiking neural signal, comprising the following steps: 1) pretreatment of spiking neural signal, 2) modeling of posterior cortical spiking neural signal generation probability, 3) model accuracy measurement, 4) model optimization with the help of numerical gradient descent and 5) prediction of posterior cortical spiking neural signal with iterative calculation. The prediction method aims to incorporate natural features of point process of spiking neural signal into the target for optimization of prediction model to improve model's capability in predicting neural spike trains.

Claims

1. A method for prediction of encephalic pulse neutral signal of a subject, characterized in that the method comprises the following steps in the following order: 1) collecting original encephalic pulse neural signal from the subject by using micro electrode array; dividing the original encephalic pulse neutral signal into a plurality of time slots by using 10 milliseconds as a time slot breadth; designating the time slot with existence of the original pulse neural signal in the time slot as 1 and designating the time slot without existence of the original pulse neural signal in the time slot as 0 to complete discretization; dividing each pulse neural signal channel into 50% training set, 20% verification set and 30% test set respectively along a time axis; checking whether pulse release rate in the training set, the verification set and the test set of the signal channel is between 2 Hz and 40 Hz; screening out signal channels with pulse release rate that is not between 2 Hz and 40 Hz; 2) using mutual information to obtain numerous anterior encephalic signal channels with the maximum value of mutual information with each signal channel in posterior encephalic region as an independent variable of a model of posterior encephalic pulse neutral signal generation probability; after that, setting sliding window in view of each time slot in all posterior encephalic signal channels to be predicted to select input independent variable of the model; once the input independent variable of the model is obtained, proceeding with modeling of posterior encephalic pulse neutral signal generation probability by using Nonhomogeneous Poisson point process generalized linear model; 3) using the model in Step 2) to obtain a sequence for pulse release probability of pulse neural signal to be predicted in the posterior encephalic region; using Discrete Time Rescaling Kolmogorov-Smirnov Statistics as a dimension for measurement of model effect to calculate model effect value; 4) taking the model effect measurement value obtained in Step 3) as a target function for derivation to obtain a numerical gradient; proceeding with optimized calculation of optimal parameters with Quasi-Newton Method; 5) executing Step 3) and 4) iteratively when calculating optimal parameters until accuracy of the model corresponding to the optimized parameters is up to a local limit point or upper limit of iterations; introducing final parameters into the model, and using the model to complete prediction of posterior encephalic pulse neutral signal; and 6) outputting predicted posterior encephalic pulse neutral signal through a display device.

2. The method for prediction of encephalic pulse neutral signal according to claim 1, characterized in that the method for obtainment independent variable of the model in the Step 2) is stated as follows: once two-two mutual information value for all pulse neutral signals of anterior and posterior encephalic regions is calculated, selecting first 8 anterior encephalic signal channels with the highest mutual information value with each pulse neural signal in posterior encephalic region as the independent variable of model for the posterior encephalic signal channel; mutual information value is stated as follows: $\begin{array}{cc}I\left(x_i,y_j\right)=\underset{x_t\in x_i}{.Math.}\underset{y_t\in y_i}{.Math.}p\left(x_t,y_t\right)\log \left(\frac{p\left(x_t,y_t\right)}{p\left(x_t\right)p\left(y_t\right)}\right)& \left(1\right)\end{array}$ wherein, x.sub.i is the signal channel i of anterior encephalic region; y.sub.j is the signal channel j of posterior encephalic region; p(x.sub.t, y.sub.t) is the joint probability for simultaneous occurrence of event x.sub.t and y.sub.t; p(x.sub.t) or p(y.sub.t) refers to occurrence probability of event x.sub.t or y.sub.t respectively; as original encephalic pulse neural signal has been discretized, value range of x.sub.t and y.sub.t is {0, 1}.

3. The method for prediction of encephalic pulse neutral signal according to claim 2, characterized in that the method for selection of input independent variable for the model is stated as follows: in view of each time slot in the posterior encephalic signal channel to be predicted, setting a sliding window with time span of 90 milliseconds to intercept the first 8 anterior encephalic pulse signal channels with the highest mutual information value with it as the input independent variable for the model; the end of this sliding window is in flush with current time slot.

4. The method for prediction of encephalic pulse neutral signal according to claim 3, characterized in that in the Step 2), Nonhomogeneous Poisson point process generalized linear model is used for modeling of generation probability of posterior encephalic pulse neutral signal; the specific form is stated as follows: $\begin{array}{cc}P\left(n\left(t\right)\right)=\frac{e^{-\lambda \left(t.Math.x\left(t\right)\right)\Delta }*{\left(\lambda \left(t.Math.x\left(t\right)\right)\Delta \right)}^{n\left(t\right)}}{n\left(t\right)!}& \left(2\right)\end{array}$ wherein, n(t) refers to the number of pulses released in time slot t; in view of discretization to original encephalic pulse neutral signal, value range of n(t) is {0,1}; P(n(t)) refers to the probability for release of n(t) pulses in time slot t; λ(t|x(t)) refers to density function for conditional probability for release of pulses in time slot t; x(t) refers to anterior encephalic pulse neutral signal collected by the sliding window; Δ refers to the length of time slot; then, linking up λ(t) and independent variable x(t), and using exponential function as the universalized linear model to link function; specific form is stated as follows:
λ(t)Δ=e.sup.x(t).sup.T.sup.θ  (3) wherein λ(t)=λ(t|x(t)).

5. The method for prediction of encephalic pulse neutral signal according to claim 4, characterized in that in Step 3), Discrete Time Rescaling Kolmogorov-Smirnov Statistics are taken as the dimension for measurement of model effect to calculate model effect value; its specific form is stated as follows: (I), firstly, generating a new sequence q(t) on the basis of λ(t) to eliminate impact from sampling of discretization time, which is represented by:
q(t)=−log(1−λ(t)Δ)  (4) (II) then, cutting the sequence q(t) in reference to the pulse sequence produced by the neural cell to be predicted; once completed, integrating value q(t) between two adjacent pulse signals, and correcting discrepancy incurred by discretization time treatment, which is represented by: $\begin{array}{cc}\zeta \left(i\right)=\underset{t=t_i+1}{\overset{t_{i+1}}{.Math.}}q\left(t\right)+q\left(t_{i+1}\right)\frac{\delta \left(i\right)}{\Delta }& \left(5\right)\end{array}$ wherein, t.sub.i refers to the occurrence time of pulse signal i of the neural cell to be predicted; t.sub.i+1 refers to occurrence time of pulse signal i+1 of the neural cell to be predicted; δ(i) is a random variable collected in the following form: $\begin{array}{cc}\delta \left(i\right)=-\frac{\Delta }{q\left(t_{i+1}\right)}\log \left(1-r\left(i\right)\left(1-e^{-q\left(t_{i+1}\right)}\right)\right)& \left(6\right)\end{array}$ r(i) in Formula (6) is a random variable uniformly distributed within [0,1]; (III) then, proceeding with scale reduction of ζ(i) to obtain the following formula:
z(i)=1−e.sup.−ζ(i)  (7) according to Kolmogorov-Smirnov Test Theory, the set {z(i)} produced by Formula (7) is to be in uniform distribution within [0,1]; (IV) rearranging the set {z(i)} in a sequence from small value to big one, and plotting it to the coordinate system for comparison while ensuring uniform distribution within standard [0,1]; transversing axis refers to set {z(i)}; whereas vertical axis refers to standard uniform distribution; taking the maximum value difference between the two on the vertical axis as the effect measurement value f(y|θ).

6. The method for prediction of encephalic pulse neutral signal according to claim 5, characterized in that in Step 3), calculation process for each neural cell signal channel to be predicted is executed for 10 times for taking the mean value when calculating effect measurement value.

7. The method for prediction of encephalic pulse neutral signal according to claim 5, characterized in that in Step 4) model optimization refers to the following fact: directly calculating numerical gradient of target function f(y|θ) through approximate derivation; the specific form for derivation of numerical gradient is stated as follows: $\begin{array}{cc}\frac{\partial f\left(y.Math.\theta \right)}{\partial {\theta }_i}=\frac{f\left(y.Math.{\theta }^a\right)-f\left(y.Math.{\theta }^b\right)}{2\mathrm{.Math.}}& \left(8\right)\end{array}$ wherein, y refers to pulse neural signal channel to be predicted in posterior encephalic region; θ refers to the vector of model parameter to be defined; f(y|θ) refers to the target function; θ.sup.a and θ.sup.b are represented as θ.sub.i.sup.a=θ.sub.i+ϵ and θ.sub.i.sup.b=θ.sub.i+ϵ on the dimension i; value of other dimensions is consistent with θ; here, value of ϵ is 1e-5.

Description

BRIEF DESCRIPTION OF DRAWINGS

(1) FIG. 1 is the overall flow chart for prediction of posterior cortical spiking neural signal in embodiments; and

(2) FIG. 2 is the graph showing dependence between each variable predicted by posterior spiking neural signal in embodiments.

PREFERRED EMBODIMENTS

(3) The present invention is further described as follows in combination with preferred embodiments and drawings.

(4) The method for prediction of cortical spiking neural signal as shown in FIGS. 1 and 2, comprising the following steps:

(5) (1) Pretreatment of Cortical Spiking Neural Signal Channel

(6) First, using micro electrode array to collect original cortical spiking neural signal, and then using time slot to divide original cortical spiking neural signal to complete discretization. Specific method is stated as follows: taking 10 milliseconds as the time slot width to divide original cortical spiking neural signal; recording the time slot with existence of spiking neural signal in the time slot as 1 or 0 on the contrary to complete discretization.

(7) After that, screening out signal channels with extremely high or low spike release rate among all cortical spiking neural signal channels; the specific method is stated as follows: First, dividing each spiking neural signal channel into 50% training set, 20% validation set and 30% test set respectively as per time axis; checking if spiking release rate in three sets of the signal channel is between 2 Hz and 40 Hz; it will be deemed to be qualified through screening if this condition can be satisfied.

(8) (2) Modeling of Posterior Cortical Spiking Neural Signal Generation Probability

(9) Once mutual information value for all pairs of spiking neural signals of anterior and posterior cortical regions is calculated, selecting the first 8 anterior cortical signal channels with the highest mutual information value with each spiking neural signal in posterior cortical region as the independent variable of model for the posterior cortical signal channel; mutual information value is stated as follows:

(10) $\begin{array}{cc}I\left(x_i,y_j\right)=\underset{x_t\in x_i}{.Math.}\underset{y_t\in y_i}{.Math.}p\left(x_t,y_t\right)\log \left(\frac{p\left(x_t,y_t\right)}{p\left(x_t\right)p\left(y_t\right)}\right)& \left(1\right)\end{array}$

(11) wherein, x.sub.i is the signal channel i of anterior cortical region; y.sub.j is the signal channel j of posterior cortical region; p(x.sub.t, y.sub.t) is the joint probability for simultaneous occurrence of event x.sub.t and y.sub.t; p(x.sub.t) or p(y.sub.t) refers to occurrence probability of event x.sub.t or y.sub.t respectively; as original cortical spiking neural signal has been discretized, value range of x.sub.t and y.sub.t is {0, 1}.

(12) In view of each time slot in the posterior cortical signal channel to be predicted, setting a sliding window with time span of 90 milliseconds to intercept the first 8 anterior cortical spiking signal channels with the highest mutual information value with it as the input independent variable for the model; the end of this sliding window is in flush with current time slot.

(13) Once the input independent variable for the model is obtained, using inhomogeneous Poisson point process generalized linear model for modeling of generation probability of posterior cortical spiking neural signal; the specific form is stated as follows:

(14) $\begin{array}{cc}P\left(n\left(t\right)\right)=\frac{e^{-\lambda \left(t.Math.x\left(t\right)\right)\Delta }*{\left(\lambda \left(t.Math.x\left(t\right)\right)\Delta \right)}^{n\left(t\right)}}{n\left(t\right)!}& \left(2\right)\end{array}$

(15) wherein, n(t) refers to the number of spikes released in time slot t; in view of discretization to original cortical spiking neural signal, value range of n(t) is {0,1}; P (n(t)) refers to the probability for release of n(t) spikes in time slot t; λ(t|x(t)) refers to density function(abbreviated as λ(t)) for conditional probability for release of spikes in time slot t; x(t) refers to anterior cortical spiking neural signal collected by the sliding window; Δ refers to the length of time slot;

(16) After that, linking up λ(t) and independent variable x(t), and use exponential function as the link function of generalized linear model; specific form is stated as follows:
λ(t)Δ=e.sup.x(t).sup.T.sup.θ  (3)

(17) wherein λ(t)=λ(t|x(t)),

(18) (3) Model Accuracy Measurement

(19) Using the model in Step 2) to obtain the series for spiking release probability of spiking neural signal to be predicted in posterior cortical region; taking Discrete Time Rescaling Kolmogorov-Smirnov Statistics as the dimension for measurement of model effect to calculate model effect measurement value; the specific form is stated as follows:

(20) I. Firstly, generating a new sequence q(t) on the basis of λ(t) to eliminate impact from sampling of discretization time; its specific form is stated as follows:
q(t)=−log(1−λ(t)Δ)  (4)

(21) II. After that, cutting the sequence q(t) in reference to the spiking trains produced by the neural cell to be predicted; once completed, integrating value q(t) between two adjacent spiking signals, and correcting discrepancy incurred by discretization time treatment; the specific form is stated as follows:

(22) $\begin{array}{cc}\zeta \left(i\right)=\underset{t=t_i+1}{\overset{t_{i+1}}{.Math.}}q\left(t\right)+q\left(t_{i+1}\right)\frac{\delta \left(i\right)}{\Delta }& \left(5\right)\end{array}$

(23) wherein, t.sub.i refers to the occurrence time of spiking signal i of the neural cell to be predicted; t.sub.i+1 refers to occurrence time of spiking signal i+1 of the neural cell to be predicted; δ(i) is a random variable collected in the following form:

(24) $\begin{array}{cc}\delta \left(i\right)=-\frac{\Delta }{q\left(t_{i+1}\right)}\log \left(1-r\left(i\right)\left(1-e^{-q\left(t_{i+1}\right)}\right)\right)& \left(6\right)\end{array}$

(25) r(i) in Formula (6) is a random variable uniformly distributed within [0,1];

(26) III. After that, proceeding with scale transformation of ζ(i) to obtain the following formula:
z(i)=1−e.sup.−ζ(i)  (7)

(27) According to Kolmogorov-Smirnov Test Theory, the set {z(i)} produced by Formula (7) is to be in uniform distribution within [0,1];

(28) IV. Rearranging the set {z(i)} in a sequence from small value to big one, and plotting it to the coordinate system for comparison while ensuring uniform distribution within standard [0,1]; transversing axis refers to set {z(i)}; whereas vertical axis refers to standard uniform distribution; taking the maximum value difference between the two on the vertical axis as the effect measurement value f(y|θ). As effect measurement value f(y|θ) may change accompanied by variation to model parameter θ, it can be deemed as a function of θ, namely target function of the model (here, y refers to a posterior cortical spiking neural signal channel to be predicted). As random variable is introduced during calculation, calculation of effect measurement value is to be executed for numerous times so as to take the mean value to eliminate impact. For each neuron signal channel to be predicted, it is necessary to execute calculation for 10 times, and take the mean value.

(29) (4) Model Optimization with the Help of Numerical Gradient Decent

(30) Taking the effect measurement value f(y|θ) of model as obtained in Step 3) as the target function to further obtain the numerical gradient through derivation of target function; proceeding with optimized calculation of optimal parameters with Quasi-Newton Method;

(31) Directly calculating numerical gradient of target function f(y|θ) through approximate derivation; the specific form for derivation of numerical gradient is stated as follows:

(32) 0$\begin{array}{cc}\frac{\partial f\left(y.Math.\theta \right)}{\partial {\theta }_i}=\frac{f\left(y.Math.{\theta }^a\right)-f\left(y.Math.{\theta }^b\right)}{2\mathrm{.Math.}}& \left(8\right)\end{array}$

(33) wherein, y refers to spiking neural signal channel to be predicted in posterior cortical region; θ refers to the vector of model parameter to be defined; f(y|θ) refers to the target function; θ.sup.a and θ.sup.b are represented as θ.sub.i.sup.a=θ.sub.i+ϵ and θ.sub.i.sup.b=θ.sub.i−ϵ on the dimension i; value of other dimensions is consistent with θ; here, value of ϵ is 1e-5.

(34) Once derivative value of the target function is obtained through approximate calculation based on numerical gradient, using Quasi-Newton Method for optimal calculation. As calculation of numerical gradient is time consuming, optimization with Quasi-Newton Method is available for quicker convergence to local optima as compared with steepest descend method. Therefore, it can minimize iteration and calculation frequency.

(35) (5) Iterative Calculation

(36) Step 3) and 4) were executed iteratively when calculating the optimal parameters until accuracy of the model corresponding to the optimized parameters is up to the local optima or upper limit of iterations. As a criteria judgment of local optima reached, target function of the model will no longer reduce accompanied by new around of iteration. The upper limit for the preset iteration frequency is 400. The parameter optimization process was directly exited without proceeding with new round of parameter iteration if upper limit for iteration frequency is reached. Final parameters were introduced into the model, and this model was used to complete prediction of posterior cortical spiking neural signal.

(37) Effect Comparison

(38) Test results of the present invention were verified by neuron spiking signal collected by Qiushi Academy for Advanced Studies of Zhejiang University from PMd and M1 cortical regions of macaque (No. B04). The signal was collected under the background that macaque was executing center-out task in four direction, of which valid length was up to 636.4 s. It comprised 127 PMd cortical neuron spiking signal channels and 102 M1 cortical neuron spiking signal channels. 96 valid PMd signal channels and 57 valid M1 signal channels were obtained through screening.

(39) The experiment effect of the present invention was compared with effect of several other spiking signal models, including linear regression model (LR), Spike-Triggered Average(STA), Spike-Triggered Convaiance (STC), conventional generalized linear model (GLM) and so on. Table 1 showed parallel comparison between this method and several other methods in terms of spiking signal prediction effect; wherein, DTR-KS is abbreviation of Discrete Time Rescaling Kolmogorov Smirnov Test, of which value is within [0,1]; the value that is more approximate to zero indicates better prediction effect of the model.

(40) TABLE-US-00001 TABLE 1 Effect Comparison between This Method and Several Other Prediction Methods M1 cortical signal channel Prediction Mean value of with the optimal effect method DTR-KS obtained with this method Linear regression 0.1395 (±0.0627) 10 model STA 0.1241 (±0.0546) 12 STC 0.1387 (±0.0865) 2 GLM 0.1176 (±0.0481) 6 This method 0.1134 (±0.0463) 27