A 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.

Description

FIELD OF THE INVENTION

The present invention is related to the spiking trains prediction filed, in particular to a method for prediction of cortical spiking trains.

BACKGROUND ARTS

Accompanied by sustainable development of technologies on micro electrode array and profound study of brain functions in the field of modern medicine, neural spiking signal can be collected and analyzed in an extremely high temporal and spatial resolution. This can facilitate scientific research personnel to explore brain signal transmission mechanism and functional link of each cerebral cortices.

The brain can be roughly divided into numerous cortical regions according to division of different functions. It is an approach to inquire into functional links between two linked cortical regions in terms of physical distribution by using cortical spiking signal from anterior cortical region to predict cortical spiking signal from posterior cortical region. Due to the working principle on which neuron takes spiking signal as the basic unit for neuron communication, cortical signal produced by single neuron can be deemed as a time series point process. Nevertheless, most of existing models have failed to incorporate point process features of cortical signal into the model for consideration. This has brought forth certain challenges to accurate prediction of spiking signal from target cortical region and judgment of functional links between two adjacent cortical regions.

SUMMARY OF THE INVENTION

In view of disadvantages to the prior arts, the present invention aims to provide a method for prediction of cortical spiking neural signal and a method for prediction of spiking neural signal that takes generalized linear model as the basis, Discrete Time Rescaling Kolmogorov-Smirnov Statistics as the optimization target and numerical gradient as optimization approach. The main problem to be solved is to incorporate natural features for point process of spiking neural signal into the optimization target for predicted model so as to improve model's capability in predicting neural spiking trains.

The technical solution provided by the present invention to solve aforesaid problems is stated as follows:

A method for prediction of cortical spiking neural signal, characterized in that it comprises the following steps:

1) Pretreatment of cortical spiking neural signal channel: firstly, using time slot to divide original cortical spiking neural signal to complete discretization; screening out signal channels with extremely high or low spiking release rate among all cortical spiking neural signal channels;

2) Modeling of posterior cortical spiking neural signal generation probability: using mutual information to obtain numerous anterior cortical signal channels with the maximum value of mutual information with each signal channel in posterior cortical region as the independent variable of the model; after that, setting sliding window in view of each time slot in all posterior cortical 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 cortical spiking neural signal generation probability with the help of inhomogeneous Poisson point process generalized linear model;

3) Model accuracy measurement: using the model in Step 2) to obtain the sequence for spiking release probability of spiking neural signal to be predicted in the posterior cortical region; taking Discrete Time Rescaling Kolmogorov-Smirnov Statistics as the dimension for measurement of model effect to calculate model effect value;

4) Model optimization with the help of numerical gradient descent: taking the model effect measurement value obtained in Step 3) as the target function for derivation to obtain the numerical gradient; proceeding with optimized calculation of optimal parameters with Quasi-Newton Method;

5) Iterative calculation: executing Step 3) and 4) 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; introducing final parameters into the model, and use this model to complete prediction of posterior cortical spiking neural signal.

Discretization in the Step 1) is stated as follows: Take 10 ms as the time slot width to divide original cortical spiking neural signal; record the time slot with existence of spiking neural signal in the time slot as 1 or 0 on the contrary to complete discretization.

The screening method in the Step 1) 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 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.

The method for obtainment independent variable of the model in the Step 2) is stated as follows: 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:

$\begin{array}{cc}I\left(x_i,y_j\right)=\sum _{x_t\in x_i}\sum _{y_t\in y_i}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_i is the signal channel i of anterior cortical region; y_j is the signal channel j of posterior cortical region; p(x_t, y_t) is the joint probability for simultaneous occurrence of event x_t and y_t ; p(x_t) or p(y_t) refers to occurrence probability of event x_t or y_t respectively; as original cortical spiking neural signal has been discretized, value range of x_t and y_t is {0, 1}.

The method for selection of input independent variable for the model in the Step 2) is stated as follows: In view of each time slot in the posterior cortical signal channel to be predicted, setting a sliding window with time span of 90ms 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.

In the Step 2), inhomogeneous Poisson point process generalized linear model is used for modeling of generation probability of posterior cortical spiking neural signal; the specific form is stated as follows:

$\begin{array}{cc}P\left(n\left(t\right)\right)=\frac{e^{-\lambda \left(tx\left(t\right)\right)\Delta }*{\left(\lambda \left(t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 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 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;

After that, linking up λ(t) and independent variable x(t), and using exponential function as the link function of generalized linear model; specific form is stated as follows:

λ(t)Δ=e^x(t)Tθ  (3)

wherein λ(t)=λ(t|x(t)).

In the 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. First, 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:

λ(t)=−log(1−λ(t)Δ)   (4)

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:

$\begin{array}{cc}\zeta \left(i\right)=\sum _{t=t_i+1}^{t_{i+1}}q\left(t\right)+q\left(t_{i+1}\right)\frac{\delta \left(i\right)}{\Delta }& \left(5\right)\end{array}$

wherein, t_j refers to the occurrence time of spiking signal i of the neural cell to be predicted; t_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:

$\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. After that, proceeding with scale transformation of ζ(i) to obtain the following formula:

z(i)=1−e^−ζ(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|θ).

In the Step 3), calculation process for each neural cell signal channel to be predicted will be executed for 10 times for taking the mean value when calculating effect measurement value.

In the 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\theta \right)}{\partial {\theta }_i}=\frac{f\left(y{\theta }^a\right)-f\left(y{\theta }^b\right)}{2\varepsilon }& \left(8\right)\end{array}$

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; θ^a and θ^b are represented as θ_i^a=θ_i+ϵ and θ_i^b=θ_i−ϵ on the dimension i; value of other dimensions is consistent with θ; here, value of ϵ is 1e-5.

As compared with prior arts, the present invention has the following beneficial effect: Natural features of point process of spiking neural signal is incorporated into optimization target for the prediction model to improve model's capability in prediction of neural spiking trains.

BRIEF DESCRIPTION OF DRAWINGS

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

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

PREFERRED EMBODIMENTS

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

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

(1) Pretreatment of Cortical Spiking Neural Signal Channel

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.

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.

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

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:

$\begin{array}{cc}I\left(x_i,y_j\right)=\sum _{x_t\in x_i}\sum _{y_t\in y_i}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_i is the signal channel i of anterior cortical region; y_j is the signal channel j of posterior cortical region; p(x_t, y_t) is the joint probability for simultaneous occurrence of event x_t and y_t; p(x_t) or p(y_t) refers to occurrence probability of event x_t or y_t respectively; as original cortical spiking neural signal has been discretized, value range of x_t and y_tis {0, 1}.

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.

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:

$\begin{array}{cc}P\left(n\left(t\right)\right)=\frac{e^{-\lambda \left(tx\left(t\right)\right)\Delta }*{\left(\lambda \left(t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 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;

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^x(t)Tθ  (3)

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

(3) Model Accuracy Measurement

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:

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)

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:

$\begin{array}{cc}\zeta \left(i\right)=\sum _{t=t_i+1}^{t_{i+1}}q\left(t\right)+q\left(t_{i+1}\right)\frac{\delta \left(i\right)}{\Delta }& \left(5\right)\end{array}$

wherein, t_i refers to the occurrence time of spiking signal i of the neural cell to be predicted; t_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:

$\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. After that, proceeding with scale transformation of ζ(i) to obtain the following formula:

z(i)=1−e^−ζ(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|θ). 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.

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

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;

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\theta \right)}{\partial {\theta }_i}=\frac{f\left(y{\theta }^a\right)-f\left(y{\theta }^b\right)}{2\varepsilon }& \left(8\right)\end{array}$

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; θ^a and θ^b are represented as θ_i^a=θ_i+ϵ and θ_i^b=θ_i−ϵ on the dimension i; value of other dimensions is consistent with θ; here, value of ϵ is 1e-5.

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.

(5) Iterative Calculation

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.

Effect Comparison

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.

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.

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

Claims

1. A method for prediction of cortical spiking neural signal, characterized in that the method comprises the following steps: 1) pretreatment of cortical spiking neural signal channel: First, using time slot to divide original cortical spiking neural signal to complete discretization; screening out signal channels with extremely high or low spiking release rate among all cortical spiking neural signal channels;2) modeling of posterior cortical spiking neural signal generation probability: using mutual information to obtain numerous anterior cortical signal channels with the maximum value of mutual information with each signal channel in posterior cortical region as the independent variable of the model; after that, setting sliding window in view of each time slot in all posterior cortical 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 cortical spiking neural signal generation probability with help of inhomogeneous Poisson point process generalized linear model;3) model accuracy measurement: using the model in Step 2) to obtain the sequence for spiking release probability of spiking neural signal to be predicted in the posterior cortical region; using Discrete Time Rescaling Kolmogorov-Smirnov Statistics as the dimension for measurement of model effect to calculate model effect value;4) model optimization with the help of numerical gradient descent: taking the model effect measurement value obtained in Step 3) as the target function for derivation to obtain the numerical gradient; proceeding with optimized calculation of optimal parameters with Quasi-Newton Method;5) iterative calculation: executing Step 3) and 4) 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; introducing final parameters into the model, and using this model to complete prediction of posterior cortical spiking neural signal.

2. The method for prediction of cortical spiking neural signal according to claim 1, characterized in that discretization in the Step 1) 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.

3. The method for prediction of cortical spiking neural signal according to claim 1, characterized in that the screening method in the Step 1) 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.

4. The method for prediction of cortical spiking neural 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 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:

5. The method for prediction of cortical spiking neural signal according to claim 4, 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 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.

6. The method for prediction of cortical spiking neural signal according to claim 5, characterized in that in the Step 2), inhomogeneous Poisson point process generalized linear model is used for modeling of generation probability of posterior cortical spiking neural signal; the specific form is stated as follows:

7. The method for prediction of cortical spiking neural signal according to claim 6, 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; its specific form is stated as follows: q(t)=−log(1−λ(t)Δ)   (4)II. then, cutting the sequence q(t) in reference to the spike train 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:

8. The method for prediction of cortical spiking neural signal according to claim 7, 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.

9. The method for prediction of cortical spiking neural signal according to claim 7, 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: