WO2018094718A1 - 一种脑区脉冲神经信号的预测方法 - Google Patents

一种脑区脉冲神经信号的预测方法 Download PDF

Info

Publication number
WO2018094718A1
WO2018094718A1 PCT/CN2016/107425 CN2016107425W WO2018094718A1 WO 2018094718 A1 WO2018094718 A1 WO 2018094718A1 CN 2016107425 W CN2016107425 W CN 2016107425W WO 2018094718 A1 WO2018094718 A1 WO 2018094718A1
Authority
WO
WIPO (PCT)
Prior art keywords
model
brain region
signal
pulse
subsequent
Prior art date
Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
Ceased
Application number
PCT/CN2016/107425
Other languages
English (en)
French (fr)
Inventor
潘纲
邢东
钱存乐
王怡雯
张巧生
郝耀耀
王跃明
郑筱祥
吴朝晖
Current Assignee (The listed assignees may be inaccurate. Google has not performed a legal analysis and makes no representation or warranty as to the accuracy of the list.)
Zhejiang University ZJU
Original Assignee
Zhejiang University ZJU
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by Zhejiang University ZJU filed Critical Zhejiang University ZJU
Priority to US16/463,686 priority Critical patent/US11238991B2/en
Publication of WO2018094718A1 publication Critical patent/WO2018094718A1/zh
Anticipated expiration legal-status Critical
Ceased legal-status Critical Current

Links

Images

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16H—HEALTHCARE INFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR THE HANDLING OR PROCESSING OF MEDICAL OR HEALTHCARE DATA
    • G16H50/00—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics
    • G16H50/50—ICT specially adapted for medical diagnosis, medical simulation or medical data mining; ICT specially adapted for detecting, monitoring or modelling epidemics or pandemics for simulation or modelling of medical disorders
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/40—Detecting, measuring or recording for evaluating the nervous system
    • A61B5/4058—Detecting, measuring or recording for evaluating the nervous system for evaluating the central nervous system
    • A61B5/4064—Evaluating the brain
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7271—Specific aspects of physiological measurement analysis
    • A61B5/7275—Determining trends in physiological measurement data; Predicting development of a medical condition based on physiological measurements, e.g. determining a risk factor
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/72—Signal processing specially adapted for physiological signals or for diagnostic purposes
    • A61B5/7221—Determining signal validity, reliability or quality
    • G—PHYSICS
    • G06—COMPUTING OR CALCULATING; COUNTING
    • G06F—ELECTRIC DIGITAL DATA PROCESSING
    • G06F2111/00—Details relating to CAD techniques
    • G06F2111/10—Numerical modelling

Definitions

  • the invention relates to the field of pulse neural signal prediction, and in particular to a method for predicting a pulsed neural signal in a brain region.
  • brain neuron pulse signals can be collected and analyzed with extremely high temporal and spatial resolution, which is for researchers to explore brain signaling mechanisms and Functional connections in the brain area are facilitated.
  • the brain can be roughly divided into several brain regions according to different functional partitions.
  • Using the brain signal of the signal brain region to predict the brain signal in the subsequent brain region is a major way to understand the functional connection between two brain regions connected in physical distribution. Since the brain neurons work with the pulse signal as the basic unit for inter-neuron communication, the brain signal generated by a single neuron can be regarded as a time series point process.
  • most of the existing models do not incorporate the point process characteristics of brain signals into the model considerations, which poses certain challenges for accurately predicting the pulse signals of the target brain regions and judging the functional connections between adjacent brain regions.
  • the object of the present invention is to provide a method for predicting a pulsed neural signal in a brain region according to the deficiencies of the prior art. Based on a general linear model, the statistic of the Kolmonov-Smirnov is reproduced in discrete time.
  • the optimization objective and the pulse neural signal prediction method with numerical gradient as the optimization method mainly solve the problem that the natural characteristics of the point process of the pulse neural signal are included in the optimization model of the prediction model. Targeting, thereby improving the predictive power of the model for neural pulse sequences.
  • a method for predicting a pulsed neural signal in a brain region includes the following steps:
  • pulsed neural signal channel preprocessing in the brain region firstly, the pulse signal of the original brain region is discretized by time slot division; the signal channel with too high or too low pulse rate in all pulse brain signal channels is screened out;
  • Model accuracy measure The pulse release probability sequence of the pulsed neural signal to be predicted in the subsequent brain region is obtained by the model in step 2); the Kolmonov-Smirnov statistic is used as a measure by discrete time replay The dimension of the model effect, calculating the effect metric of the model;
  • step 4) Using the numerical gradient descent to optimize the model: taking the effect metric of the model obtained in step 3) as the objective function, then deriving the objective function to obtain the numerical gradient, and using the quasi-Newton method to optimize the optimal parameters;
  • the discretization in the step 1) refers to dividing the pulse signal of the original brain region by using 10 milliseconds as the time slot width, and recording the time slot of the pulse neural signal in the time slot as 1 and vice versa. Complete the discretization.
  • the screening method in the step 1) refers to: firstly dividing each pulse neural signal channel into 50% training set, 20% verification set and 30% test set according to the time axis, and detecting the pulse in the three sets of the signal channel. Whether the issuance rate is between 2Hz and 40Hz, and the screening is met if this condition is met.
  • the highest first 8 channels of the brain signal channel are used as model independent variables of the signal channel of the subsequent brain region.
  • the mutual information value is:
  • x i is the i-th signal path of the anterior brain region
  • y j is the j-th signal path of the succeeding brain region
  • p(x t , y t ) is the joint probability of the simultaneous occurrence of the event x t and y t
  • p( x t ) and p(y t ) represent the probability of occurrence of the events x t and y t respectively. Since the original brain region pulsed neural signals have been discretized, the values of x t and y t are in the range ⁇ 0, 1 ⁇ . between.
  • the general linear model of the non-homogeneous Poisson point process is used to model the probability model of the pulsed neural signal generation in the subsequent brain region, and the specific form is as follows:
  • n(t) represents the number of pulses issued in the t-th time slot. Since the pulse signal of the original brain region is discretized, the value of n(t) is ⁇ 0, 1 ⁇ ; P(n( t)) indicates the probability of issuing n(t) pulses in the tth time slot; ⁇ (t
  • ⁇ (t) is associated with the independent variable x(t), and a general linear model with an exponential function as a connection function is adopted.
  • the specific form is:
  • the discrete Körmonov-Smirnov statistic is used as the dimension of the measurement model effect, and the effect metric value of the model is calculated, and the specific form is:
  • sequence q(t) is cut according to the pulse sequence generated by the neuron to be predicted, and then the q(t) value between the adjacent two pulse signals is integrated, and the deviation generated by the discrete time processing is performed. Corrected, the specific form is:
  • t i represents the time at which the ith pulse signal of the predicted neuron occurs
  • t i+1 refers to the time at which the i+1th pulse signal of the neuron to be predicted occurs
  • ⁇ (i) is acquired in the following form The resulting random variable:
  • r(i) in equation (6) is a random variable uniformly distributed between [0, 1];
  • the set ⁇ z(i) ⁇ is rearranged according to the numerical value from small to large, and is uniformly plotted with the standard [0, 1] in the coordinate system for comparison, and the horizontal axis is the set ⁇ z(i) ) ⁇ , the vertical axis is a standard uniform distribution, and the maximum numerical difference between the two is taken as the effect metric f(y
  • the calculation process is performed 10 times and averaged for each signal path to be predicted.
  • the model optimization in the step 4) refers to directly calculating the numerical gradient of the objective function f(y
  • ⁇ is the model parameter vector to be determined
  • ⁇ ) is the objective function
  • ⁇ a and ⁇ b are respectively expressed in the i-th dimension as as well as The values of the remaining dimensions are consistent with ⁇ , where ⁇ is 1e-5.
  • the beneficial effects of the present invention are embodied in: incorporating the natural characteristics of the point process of the pulsed neural signal into the optimization target of the prediction model, thereby improving the prediction ability of the model for the neural pulse sequence.
  • 1 is an overall flow chart of pulsed neural signal prediction in a subsequent brain region in an embodiment
  • Figure 2 is a graph of dependence between variables in the prediction of subsequent brain region pulsed neural signals in the examples.
  • the method for predicting pulsed neural signals in the brain region is as follows:
  • the micro-electrode array is used to collect the pulse signal of the original brain region, and then the pulse signal of the original brain region is discretized by time slot division.
  • the specific method 10 milliseconds as the time slot width is used to equidistantly divide the pulse signal of the original brain region.
  • the time slot in which the pulse neural signal exists in the time slot is recorded as 1, and the negative is recorded as 0, thereby completing the discretization.
  • each pulse neural signal channel is divided into 50% training set, 20% verification set and 30 according to the time axis.
  • the % test set detects whether the burst rates in the three sets of this signal channel are between 2 Hz and 40 Hz. After this condition is met, the filter is passed.
  • the signal path of the first 8 channels of the anterior brain region with the highest mutual information value is selected as the pulse signal of each pulse in the subsequent brain region.
  • the mutual information value is:
  • x i is the i-th signal path of the anterior brain region
  • y j is the j-th signal path of the succeeding brain region
  • p(x t , y t ) is the joint probability of the simultaneous occurrence of the event x t and y t
  • p( x t ) and p(y t ) represent the probability of occurrence of the events x t and y t respectively. Since the original brain region pulsed neural signals have been discretized, the values of x t and y t are in the range ⁇ 0, 1 ⁇ . between.
  • a sliding window with a time span of 90 milliseconds is set for each time slot in each subsequent brain region signal channel to be predicted, and the first 8 channels of the preceding brain region pulse signal channel with the highest mutual information value are intercepted.
  • the end of this sliding window is flush with the current time slot.
  • n(t) represents the number of pulses issued in the t-th time slot. Since the pulse signal of the original brain region is discretized, the value of n(t) is ⁇ 0, 1 ⁇ ; P(n( t)) indicates the probability of issuing n(t) pulses in the tth time slot; ⁇ (t
  • ⁇ (t) is associated with the independent variable x(t), and a general linear model with an exponential function as a connection function is adopted.
  • the specific form is:
  • the pulse release probability sequence of the pulsed neural signal to be predicted in the subsequent brain region is obtained by the model in step 2); the discrete Körmonov-Smirnov statistic is used as the dimension of the measurement model effect, and the calculation model is calculated by using discrete time.
  • the effect metric the specific form is:
  • sequence q(t) is cut according to the pulse sequence generated by the neuron to be predicted, and then the q(t) value between the adjacent two pulse signals is integrated, and the deviation generated by the discrete time processing is performed. Corrected, the specific form is:
  • t i represents the time at which the ith pulse signal of the predicted neuron occurs
  • t i+1 refers to the time at which the i+1th pulse signal of the neuron to be predicted occurs
  • ⁇ (i) is acquired in the following form The resulting random variable:
  • r(i) in equation (6) is a random variable uniformly distributed between [0, 1];
  • the set ⁇ z(i) ⁇ is rearranged according to the numerical value from small to large, and is uniformly plotted with the standard [0, 1] in the coordinate system for comparison, and the horizontal axis is the set ⁇ z(i) ) ⁇ , the vertical axis is a standard uniform distribution, and the maximum numerical difference between the two is taken as the effect metric f(y
  • ⁇ ) varies with the parameter ⁇ of the model, so it can be regarded as a function of ⁇ , that is, the objective function of the model (where y is the pulsed neural signal of the brain region to be predicted) aisle). Since random variables are introduced in the calculation process, the calculation process of this effect metric needs to be performed several times and averaged to eliminate the influence. For each neuron signal path to be predicted, the calculation process is performed 10 times and averaged.
  • the effect metric of the model obtained in step 3) is f(y
  • ⁇ ) is directly calculated by approximating the objective function f(y
  • the specific form of the numerical gradient is:
  • ⁇ is the model parameter vector to be determined
  • ⁇ ) is the objective function
  • ⁇ a and ⁇ b are respectively expressed in the i-th dimension as as well as The values of the remaining dimensions are consistent with ⁇ , where ⁇ is 1e-5.
  • the quasi-Newton method is used to optimize the calculation. Since the numerical gradient calculation process is time consuming, the optimization using the quasi-Newton method can converge to the local extremum faster than the steepest descent method, thereby minimizing the number of iterations, so as to achieve the purpose of reducing the calculation amount.
  • step 3) and step 4) are iteratively performed until the accuracy of the model corresponding to the optimized parameter reaches the local extreme point or reaches the upper limit of the number of iterations.
  • the criterion for reaching the local extreme point is that the objective function value of the model is no longer reduced with a new iteration.
  • the upper limit of the preset number of iterations is 400. If the upper limit of the number of iterations has been reached, the parameter optimization process is directly exited. Perform a new round of parameter iteration. The final parameters are substituted into the model, and the model is used to predict the pulsed neural signals in the subsequent brain regions.
  • test results of the present invention were tested on the neuron pulse signals of the PMd and M1 brain regions of the macaques (numbered B04) collected by the Institute of Advanced Studies of Zhejiang University.
  • the background of the signal is that the macaque performs the four-direction center-out.
  • the effective time of the signal is 636.4 seconds, including the neuron pulse signal path of the 127-channel PMd brain region and the 102-channel M1 brain region. After screening, 96 effective PMd signal pathways and 57 effective M1 signaling pathways are obtained.
  • the experimental results of the present invention are simultaneously compared with those obtained by several other pulse signal models, including linear regression model (LR), Spike-Triggered Average (STA), and Spike-Triggered. Convaiance (STC), traditional general linear model (GLM), etc.
  • Table 1 is a horizontal comparison of the prediction effects of the method and the remaining types of methods for pulse signals, where DTR-KS is the abbreviation for Discrete Time Rescaling Kolmogorov Smirnov test. The value is between [0, 1], and the closer the value is to 0, the better the prediction effect of the model.

Landscapes

  • Health & Medical Sciences (AREA)
  • Engineering & Computer Science (AREA)
  • Medical Informatics (AREA)
  • Public Health (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • General Health & Medical Sciences (AREA)
  • Pathology (AREA)
  • Biomedical Technology (AREA)
  • Databases & Information Systems (AREA)
  • Epidemiology (AREA)
  • Data Mining & Analysis (AREA)
  • Primary Health Care (AREA)
  • Neurology (AREA)
  • Animal Behavior & Ethology (AREA)
  • Veterinary Medicine (AREA)
  • Physics & Mathematics (AREA)
  • Surgery (AREA)
  • Biophysics (AREA)
  • Heart & Thoracic Surgery (AREA)
  • Molecular Biology (AREA)
  • Physiology (AREA)
  • Artificial Intelligence (AREA)
  • Computer Vision & Pattern Recognition (AREA)
  • Psychiatry (AREA)
  • Signal Processing (AREA)
  • Neurosurgery (AREA)
  • Psychology (AREA)
  • Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
  • Magnetic Resonance Imaging Apparatus (AREA)
  • Measuring And Recording Apparatus For Diagnosis (AREA)
  • Image Analysis (AREA)

Abstract

一种脑区脉冲神经信号的预测方法，包括如下步骤：1)、脑区脉冲神经信号通道预处理；2)、后继脑区脉冲神经信号生成概率模型建模；3)、模型准确性度量；4)、利用数值梯度下降进行模型优化；5)迭代计算完成后继脑区脉冲神经信号预测。本预测方法将脉冲神经信号的点过程的天然特性纳入预测模型的优化目标中，从而提高模型对神经脉冲序列的预测能力。

Description

一种脑区脉冲神经信号的预测方法 技术领域
本发明涉及脉冲神经信号预测领域，具体涉及一种脑区脉冲神经信号的预测方法。
背景技术
随着微电极阵列技术的不断发展和现代医学对脑功能研究的不断深入，大脑神经元脉冲信号能够以极高的时空分辨率被采集和分析，这为科研人员们探索大脑信号传递机制和各个脑区的功能联系提供了便利。
大脑可以依照不同的功能分区粗略划分为若干脑区。利用信号前继脑区的脑信号来预测信号后继脑区的脑信号，是了解两个在物理分布上相连的脑区之间功能联系的一大途径。由于大脑神经元以脉冲信号为基本单位进行神经元间通信的工作原理，单个神经元所产生的脑信号本身可视为一个时间序列点过程。然而，现有的多数模型大多未将脑信号的点过程特性纳入模型的考量中，这对准确预测目标脑区的脉冲信号、判断相邻脑区之间的功能联系带来一定的挑战。
发明内容
本发明的目的在于针对现有技术的不足，提供一种脑区脉冲神经信号的预测方法，以通用线性模型为基础、以离散时间重放缩柯尔莫诺夫-斯米尔诺夫统计量为优化目标、以数值梯度为优化手段的脉冲神经信号预测方法，主要解决的问题是将脉冲神经信号的点过程的天然特性纳入预测模型的优化目 标中，从而提高模型对神经脉冲序列的预测能力。
本发明解决上述技术问题所提供的技术方案为：
一种脑区脉冲神经信号的预测方法，包括如下步骤：
1)、脑区脉冲神经信号通道预处理：首先将原始脑区脉冲神经信号采用时间槽划分完成离散化；筛除所有脑区脉冲神经信号通道中脉冲发放率过高或过低的信号通道；
2)、后继脑区脉冲神经信号生成概率模型建模：利用互信息获取与后继脑区的每一路信号通道互信息值最高的若干前继脑区信号通道作为模型自变量，然后针对每一个待预测的后继脑区信号通道中的每一个时间槽设置滑动窗口来选取模型的输入自变量；得到模型的输入自变量之后，以非齐次泊松点过程通用线性模型进行后继脑区脉冲神经信号生成概率模型建模；
3)、模型准确性度量：通过步骤2)中的模型得到后继脑区待预测脉冲神经信号的脉冲发放概率序列；利用离散时间重放缩柯尔莫诺夫-斯米尔诺夫统计量作为度量模型效果的量纲，计算模型的效果度量值；
4)、利用数值梯度下降进行模型优化：将步骤3)中得到的模型的效果度量值作为目标函数，然后对目标函数求导得到数值梯度，采用拟牛顿法进行最优参数优化计算；
5)迭代计算：在求解最优参数计算的过程中，迭代执行步骤3)和步骤4)，直至优化后的参数对应的模型准确性达到局部极值点或达到迭代次数上限，将最终确定的参数代入模型中，并利用此模型完成后继脑区脉冲神经信号预测。
所述步骤1)中离散化是指：采取10毫秒作为时间槽宽度对原始脑区脉冲神经信号进行划分，将时间槽内存在脉冲神经信号的时间槽记为1，反之记为0，以此完成离散化。
所述步骤1)中筛除方法是指：首先将每一路脉冲神经信号通道按照时间轴划分成50％训练集、20％验证集和30％测试集，检测此信号通道的三个集合中脉冲发放率是否均在2Hz至40Hz之间，满足此条件便通过筛选。
所述步骤2)中获取模型自变量的方法：计算前继脑区与后继脑区各路脉冲神经信号两两互信息值之后，针对后继脑区每一路脉冲神经信号，选取与之互信息值最高的前8路前继脑区信号通道作为此后继脑区信号通道的模型自变量。互信息值为：
Figure PCTCN2016107425-appb-000001
其中，xi为前继脑区第i路信号通路，yj为后继脑区第j路信号通路，p(xt,yt)是事件xt与yt同时发生的联合概率，p(xt)和p(yt)分别表示事件xt与yt各自发生的概率，由于原始脑区脉冲神经信号已经被离散化，因此xt和yt的取值范围在{0，1}之间。
所述步骤2)中选取模型的输入自变量的方法：针对每一个待预测的后继脑区信号通道中的每一个时间槽设置一个时间跨度为90毫秒的滑动窗口，用来截取与之互信息值最高的前8路前继脑区脉冲信号通道，作为模型的输入自变量，此滑动窗口的末端与当前时间槽保持齐平。
所述步骤2)中以非齐次泊松点过程通用线性模型进行后继脑区脉冲神经信号生成概率模型建模，其具体形式为：
Figure PCTCN2016107425-appb-000002
其中，n(t)表示第t个时间槽内发放脉冲的个数，由于原始脑区脉冲神经信号采取离散化处理，因此n(t)取值范围为{0,1}；P(n(t))表示在第t个时间槽内发放n(t)个脉冲的概率；λ(t|x(t))表示第t个时间槽内发放脉冲的条件概率密 度函数；x(t)是由滑动窗口采集得到的前继脑区脉冲神经信号，Δ是时间槽的长度；
随后将λ(t)与自变量x(t)联系在一起，采用以指数函数为连接函数的通用线性模型，具体形式为：
Figure PCTCN2016107425-appb-000003
其中，λ(t)＝λ(t|x(t))。
所述步骤3)中利用离散时间重放缩柯尔莫诺夫-斯米尔诺夫统计量作为度量模型效果的量纲，计算模型的效果度量值，其具体形式是：
I、首先为了消除离散时间采样带来的影响，在λ(t)的基础之上生成一个新序列q(t)，其具体形式为：
q(t)＝-log(1-λ(t)Δ)       (4)
II、随后将序列q(t)按照待预测神经元产生的脉冲序列进行切割，然后对相邻的两个脉冲信号之间的q(t)值进行积分，并对离散时间处理产生的偏差进行修正，具体形式为：
Figure PCTCN2016107425-appb-000004
其中，ti指待预测神经元的第i个脉冲信号发生的时间，ti+1指待预测神经元的第i+1个脉冲信号发生的时间，δ(i)是一个按照以下形式采集得到的随机变量：
Figure PCTCN2016107425-appb-000005
式(6)中的r(i)是一个在[0,1]之间均匀分布的随机变量；
III、随后对ζ(i)进行尺度放缩，得到：
z(i)＝1-e-ζ(i)      (7)
根据柯尔莫诺夫-斯米尔诺夫检验理论，式(7)产生的集合{z(i)}的分布应 遵循在[0,1]之间的均匀分布；
IV、将集合{z(i)}按照数值由小到大进行重排列，并与标准[0,1]之间的均匀分布同时绘图于坐标系中进行比较，横轴为集合{z(i)}，纵轴为标准均匀分布，取纵轴上二者最大数值差作为效果度量值f(y|θ)。
所述步骤3)中计算模型的效果度量值时，针对每一个待预测神经元信号通路，计算过程均执行10遍并取平均值。
所述步骤4)中模型优化是指：通过对目标函数f(y|θ)进行近似求导直接计算目标函数f(y|θ)的数值梯度，求解数值梯度的具体形式为：
Figure PCTCN2016107425-appb-000006
其中，y是后继脑区待预测脉冲神经信号通道，θ是待确定的模型参数向量，f(y|θ)是目标函数，θa与θb在第i维上分别表示为
Figure PCTCN2016107425-appb-000007
以及
Figure PCTCN2016107425-appb-000008
其余维度的值与θ保持一致，此处ε的取值为1e-5。
同现有技术相比，本发明的有益效果体现在：将脉冲神经信号的点过程的天然特性纳入预测模型的优化目标中，从而提高模型对神经脉冲序列的预测能力。
附图说明
图1是实施例中后继脑区脉冲神经信号预测的整体流程图；
图2是实施例中后继脑区脉冲神经信号预测的各个变量之间的依赖图。
具体实施方式
下面结合具体的实施例以及附图对本发明作进一步的说明。
如图1和2所示，脑区脉冲神经信号的预测方法，具体步骤如下：
(1)脑区脉冲神经信号通道预处理
首先使用微电极阵列采集原始脑区脉冲神经信号，然后将原始脑区脉冲神经信号采用时间槽划分完成离散化，具体方法：采取10毫秒作为时间槽宽度对原始脑区脉冲神经信号进行等距划分，将时间槽内存在脉冲神经信号的时间槽记为1，反之记为0，以此完成离散化。
随后筛除所有脑区脉冲神经信号通道中脉冲发放率过高或过低的信号通道；具体方法：首先将每一路脉冲神经信号通道按照时间轴划分成50％训练集、20％验证集和30％测试集，检测此信号通道的三个集合中脉冲发放率是否均在2Hz至40Hz之间，满足此条件便通过筛选。
(2)后继脑区脉冲神经信号生成概率模型建模
计算前继脑区与后继脑区各路脉冲神经信号两两互信息值之后，针对后继脑区每一路脉冲神经信号，选取与之互信息值最高的前8路前继脑区信号通道作为此后继脑区信号通道的模型自变量。互信息值为：
Figure PCTCN2016107425-appb-000009
其中，xi为前继脑区第i路信号通路，yj为后继脑区第j路信号通路，p(xt,yt)是事件xt与yt同时发生的联合概率，p(xt)和p(yt)分别表示事件xt与yt各自发生的概率，由于原始脑区脉冲神经信号已经被离散化，因此xt和yt的取值范围在{0，1}之间。
然后针对每一个待预测的后继脑区信号通道中的每一个时间槽设置一个时间跨度为90毫秒的滑动窗口，用来截取与之互信息值最高的前8路前继脑区脉冲信号通道，作为模型的输入自变量，此滑动窗口的末端与当前时间槽保持齐平。
得到模型的输入自变量之后，以非齐次泊松点过程通用线性模型进行后 继脑区脉冲神经信号生成概率模型建模，其具体形式为：
Figure PCTCN2016107425-appb-000010
其中，n(t)表示第t个时间槽内发放脉冲的个数，由于原始脑区脉冲神经信号采取离散化处理，因此n(t)取值范围为{0,1}；P(n(t))表示在第t个时间槽内发放n(t)个脉冲的概率；λ(t|x(t))表示第t个时间槽内发放脉冲的条件概率密度函数(简写为λ(t))；x(t)是由滑动窗口采集得到的前继脑区脉冲神经信号，Δ是时间槽的长度；
随后将λ(t)与自变量x(t)联系在一起，采用以指数函数为连接函数的通用线性模型，具体形式为：
Figure PCTCN2016107425-appb-000011
其中，λ(t)＝λ(t|x(t))。
(3)模型准确性度量
通过步骤2)中的模型得到后继脑区待预测脉冲神经信号的脉冲发放概率序列；利用离散时间重放缩柯尔莫诺夫-斯米尔诺夫统计量作为度量模型效果的量纲，计算模型的效果度量值，其具体形式是：
I、首先为了消除离散时间采样带来的影响，在λ(t)的基础之上生成一个新序列q(t)，其具体形式为：
q(t)＝-log(1-λ(t)Δ)      (4)
II、随后将序列q(t)按照待预测神经元产生的脉冲序列进行切割，然后对相邻的两个脉冲信号之间的q(t)值进行积分，并对离散时间处理产生的偏差进行修正，具体形式为：
Figure PCTCN2016107425-appb-000012
其中，ti指待预测神经元的第i个脉冲信号发生的时间，ti+1指待预测神经元的第i+1个脉冲信号发生的时间，δ(i)是一个按照以下形式采集得到的随机变量：
Figure PCTCN2016107425-appb-000013
式(6)中的r(i)是一个在[0,1]之间均匀分布的随机变量；
III、随后对ζ(i)进行尺度放缩，得到：
z(i)＝1-e-ζ(i)     (7)
根据柯尔莫诺夫-斯米尔诺夫检验理论，式(7)产生的集合{z(i)}的分布应遵循在[0,1]之间的均匀分布；
IV、将集合{z(i)}按照数值由小到大进行重排列，并与标准[0,1]之间的均匀分布同时绘图于坐标系中进行比较，横轴为集合{z(i)}，纵轴为标准均匀分布，取纵轴上二者最大数值差作为效果度量值f(y|θ)。效果度量值f(y|θ)随着模型的参数θ的变化而变化，因此可看作是关于θ的一个函数，也即模型的目标函数(此处y是待预测后继脑区脉冲神经信号通道)。由于在计算过程中引入了随机变量，因此对此效果度量值的计算过程需要进行若干遍并取平均值来消除影响。针对每一个待预测神经元信号通路，计算过程均执行10遍并取平均值。
(4)利用数值梯度下降进行模型优化
将步骤3)中得到的模型的效果度量值作f(y|θ)为目标函数，然后对目标函数求导得到数值梯度，采用拟牛顿法进行最优参数优化计算；
通过对目标函数f(y|θ)进行近似求导直接计算目标函数f(y|θ)的数值梯度，求解数值梯度的具体形式为：
Figure PCTCN2016107425-appb-000014
其中，y是后继脑区待预测脉冲神经信号通道，θ是待确定的模型参数向量，f(y|θ)是目标函数，θa与θb在第i维上分别表示为
Figure PCTCN2016107425-appb-000015
以及
Figure PCTCN2016107425-appb-000016
其余维度的值与θ保持一致，此处ε的取值为1e-5。
利用数值梯度近似求得目标函数的导数值之后，采用拟牛顿法进行最优化计算。由于数值梯度计算过程较为耗时，利用拟牛顿法进行优化相比于最速下降法而言，能够更快收敛至局部极值，从而尽量减少迭代次数，以达成减少计算量的目的。
(5)迭代计算
在求解最优参数计算的过程中，迭代执行步骤3)和步骤4)，直至优化后的参数对应的模型准确性达到局部极值点或达到迭代次数上限。达到局部极值点的标准是模型的目标函数值不再随着新一轮迭代而降低，预设的迭代次数上限是400次，若已达到迭代次数上限，则直接退出参数优化过程，不再进行新一轮参数迭代。将最终确定的参数代入模型中，并利用此模型完成后继脑区脉冲神经信号预测。
效果比较
本发明的试验结果在浙江大学求是高等研究院所采集的猕猴(编号为B04)的PMd和M1脑区的神经元脉冲信号上进行了实验，该信号的背景是猕猴执行四方向center-out任务时采集而得，信号有效时长为636.4秒，包括127路PMd脑区和102路M1脑区的神经元脉冲信号通路，经过筛选后得到96路有效PMd信号通路和57路有效M1信号通路。
本发明的实验效果同时与其他几种脉冲信号模型得到的效果进行比较，包括线性回归模型(LR)、Spike-Triggered Average(STA)、Spike-Triggered  Convaiance(STC)、传统通用线性模型(GLM)等。表1是本方法与其余几类方法对于脉冲信号预测效果的横向比较，其中DTR-KS是离散时间重放缩柯尔莫诺夫-斯米尔诺夫检验(Discrete Time Rescaling Kolmogorov Smirnov test)的简称，其值在[0,1]之间，数值越接近0表示该模型的预测效果越好。
表1本方法与其他几类预测方法的效果比较
Figure PCTCN2016107425-appb-000017

Claims (9)

  1. 一种脑区脉冲神经信号的预测方法，其特征在于，包括如下步骤：
    1)、脑区脉冲神经信号通道预处理：首先将原始脑区脉冲神经信号采用时间槽划分完成离散化；筛除所有脑区脉冲神经信号通道中脉冲发放率过高或过低的信号通道；
    2)、后继脑区脉冲神经信号生成概率模型建模：利用互信息获取与后继脑区的每一路信号通道互信息值最高的若干前继脑区信号通道作为模型自变量，然后针对每一个待预测的后继脑区信号通道中的每一个时间槽设置滑动窗口来选取模型的输入自变量；得到模型的输入自变量之后，以非齐次泊松点过程通用线性模型进行后继脑区脉冲神经信号生成概率模型建模；
    3)、模型准确性度量：通过步骤2)中的模型得到后继脑区待预测脉冲神经信号的脉冲发放概率序列；利用离散时间重放缩柯尔莫诺夫-斯米尔诺夫统计量作为度量模型效果的量纲，计算模型的效果度量值；
    4)、利用数值梯度下降进行模型优化：将步骤3)中得到的模型的效果度量值作为目标函数，然后对目标函数求导得到数值梯度，采用拟牛顿法进行最优参数优化计算；
    5)迭代计算：在求解最优参数计算的过程中，迭代执行步骤3)和步骤4)，直至优化后的参数对应的模型准确性达到局部极值点或达到迭代次数上限，将最终确定的参数代入模型中，并利用此模型完成后继脑区脉冲神经信号预测。
  2. 根据权利要求1所述的脑区脉冲神经信号的预测方法，其特征在于，所述步骤1)中离散化是指：采取10毫秒作为时间槽宽度对原始脑区脉冲神经信号进行划分，将时间槽内存在脉冲神经信号的时间槽记为1，反之记为0，以此完成离散化。
  3. 根据权利要求1所述的脑区脉冲神经信号的预测方法，其特征在于，所述步骤1)中筛除方法是指：首先将每一路脉冲神经信号通道按照时间轴划分成50％训练集、20％验证集和30％测试集，检测此信号通道的三个集合中脉冲发放率是否均在2Hz至40Hz之间，满足此条件便通过筛选。
  4. 根据权利要求1所述的脑区脉冲神经信号的预测方法，其特征在于，所述步骤2)中获取模型自变量的方法：计算前继脑区与后继脑区各路脉冲神经信号两两互信息值之后，针对后继脑区每一路脉冲神经信号，选取与之互信息值最高的前8路前继脑区信号通道作为此后继脑区信号通道的模型自变量；互信息值为：
    Figure PCTCN2016107425-appb-100001
    其中，xi为前继脑区第i路信号通路，yj为后继脑区第j路信号通路，p(xt,yt)是事件xt与yt同时发生的联合概率，p(xt)和p(yt)分别表示事件xt与yt各自发生的概率，由于原始脑区脉冲神经信号已经被离散化，因此xt和yt的取值范围均在{0，1}之间。
  5. 根据权利要求4所述的脑区脉冲神经信号的预测方法，其特征在于，所述步骤2)中选取模型的输入自变量的方法：针对每一个待预测的后继脑区信号通道中的每一个时间槽设置一个时间跨度为90毫秒的滑动窗口，用来截取与之互信息值最高的前8路前继脑区脉冲信号通道，作为模型的输入自变量，此滑动窗口的末端与当前时间槽保持齐平。
  6. 根据权利要求5所述的脑区脉冲神经信号的预测方法，其特征在于，所述步骤2)中以非齐次泊松点过程通用线性模型进行后继脑区脉冲神经信号生成概率模型建模，其具体形式为：
    Figure PCTCN2016107425-appb-100002
    其中，n(t)表示第t个时间槽内发放脉冲的个数，由于原始脑区脉冲神经信号采取离散化处理，因此n(t)取值范围为{0,1}；P(n(t))表示在第t个时间槽内发放n(t)个脉冲的概率；λ(t|x(t))表示第t个时间槽内发放脉冲的条件概率密度函数；x(t)是由滑动窗口采集得到的前继脑区脉冲神经信号，Δ是时间槽的长度；
    随后将λ(t)与自变量x(t)联系在一起，采用以指数函数为连接函数的通用线性模型，具体形式为：
    Figure PCTCN2016107425-appb-100003
    其中，λ(t)＝λ(t|x(t))。
  7. 根据权利要求6所述的脑区脉冲神经信号的预测方法，其特征在于，所述步骤3)中利用离散时间重放缩柯尔莫诺夫-斯米尔诺夫统计量作为度量模型效果的量纲，计算模型的效果度量值，其具体形式是：
    I、首先为了消除离散时间采样带来的影响，在λ(t)的基础之上生成一个新序列q(t)，其具体形式为：
    q(t)＝-log(1-λ(t)Δ)   (4)
    II、随后将序列q(t)按照待预测神经元产生的脉冲序列进行切割，然后对相邻的两个脉冲信号之间的q(t)值进行积分，并对离散时间处理产生的偏差进行修正，具体形式为：
    Figure PCTCN2016107425-appb-100004
    其中，ti指待预测神经元的第i个脉冲信号发生的时间，ti+1指待预测神经元的第i+1个脉冲信号发生的时间，δ(i)是一个按照以下形式采集得到的随机 变量：
    Figure PCTCN2016107425-appb-100005
    式(6)中的r(i)是一个在[0,1]之间均匀分布的随机变量；
    III、随后对ζ(i)进行尺度放缩，得到：
    z(i)＝1-e-ζ(i)   (7)根据柯尔莫诺夫-斯米尔诺夫检验理论，式(7)产生的集合{z(i)}的分布应遵循在[0,1]之间的均匀分布；
    IV、将集合{z(i)}按照数值由小到大进行重排列，并与标准[0,1]之间的均匀分布同时绘图于坐标系中进行比较，横轴为集合{z(i)}，纵轴为标准均匀分布，取纵轴上二者最大数值差作为效果度量值f(y|θ)。
  8. 根据权利要求7所述的脑区脉冲神经信号的预测方法，其特征在于，所述步骤3)中计算模型的效果度量值时，针对每一个待预测神经元信号通路，计算过程均执行10遍并取平均值。
  9. 根据权利要求7所述的脑区脉冲神经信号的预测方法，其特征在于，所述步骤4)中模型优化是指：通过对目标函数f(y|θ)进行近似求导直接计算目标函数f(y|θ)的数值梯度，求解数值梯度的具体形式为：
    Figure PCTCN2016107425-appb-100006
    其中，y是后继脑区待预测脉冲神经信号通道，θ是待确定的模型参数向量，f(y|θ)是目标函数，θa与θb在第i维上分别表示为
    Figure PCTCN2016107425-appb-100007
    以及
    Figure PCTCN2016107425-appb-100008
    其余维度的值与θ保持一致，此处ε的取值为1e-5。
PCT/CN2016/107425 2016-11-24 2016-11-28 一种脑区脉冲神经信号的预测方法 Ceased WO2018094718A1 (zh)

Priority Applications (1)

Application Number Priority Date Filing Date Title
US16/463,686 US11238991B2 (en) 2016-11-24 2016-11-28 Method for prediction of cortical spiking trains

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
CN201611051555.4A CN106529186B (zh) 2016-11-24 2016-11-24 一种脑区脉冲神经信号的预测方法
CN201611051555.4 2016-11-24

Publications (1)

Publication Number Publication Date
WO2018094718A1 true WO2018094718A1 (zh) 2018-05-31

Family

ID=58357085

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/CN2016/107425 Ceased WO2018094718A1 (zh) 2016-11-24 2016-11-28 一种脑区脉冲神经信号的预测方法

Country Status (3)

Country Link
US (1) US11238991B2 (zh)
CN (1) CN106529186B (zh)
WO (1) WO2018094718A1 (zh)

Families Citing this family (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
WO2021072125A1 (en) * 2019-10-09 2021-04-15 University Of Southern California Geometric paradigm for nonlinear modeling and control of neural dynamics
CN111311321B (zh) * 2020-02-14 2021-11-02 北京百度网讯科技有限公司 用户消费行为预测模型训练方法、装置、设备及存储介质
US20230206038A1 (en) * 2020-05-25 2023-06-29 Zhejiang University A modeling method for artificial neural pathway across encephalic regions
CN117357134B (zh) * 2023-12-08 2024-02-09 中国科学院深圳先进技术研究院 一种神经电脉冲检测方法、系统及终端

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102360501A (zh) * 2011-07-19 2012-02-22 中国科学院自动化研究所 一种基于核磁共振成像的脑区有效连接度构建方法
CN103584855A (zh) * 2013-10-24 2014-02-19 燕山大学 脑肌电同步采集及信息传递特性分析方法
WO2015030606A2 (en) * 2013-08-26 2015-03-05 Auckland University Of Technology Improved method and system for predicting outcomes based on spatio / spectro-temporal data

Family Cites Families (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
US20030105409A1 (en) * 2001-11-14 2003-06-05 Donoghue John Philip Neurological signal decoding
CN105574324A (zh) * 2015-12-10 2016-05-11 浙江大学 一种自适应的脑神经信号处理方法及系统
US20200298005A1 (en) * 2017-05-26 2020-09-24 Newton Howard Brain-machine interface (bmi)
WO2019035887A1 (en) * 2017-08-16 2019-02-21 University Of Southern California APPARATUS AND METHOD FOR DECODING AND RESTORING COGNITIVE FUNCTIONS

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN102360501A (zh) * 2011-07-19 2012-02-22 中国科学院自动化研究所 一种基于核磁共振成像的脑区有效连接度构建方法
WO2015030606A2 (en) * 2013-08-26 2015-03-05 Auckland University Of Technology Improved method and system for predicting outcomes based on spatio / spectro-temporal data
CN103584855A (zh) * 2013-10-24 2014-02-19 燕山大学 脑肌电同步采集及信息传递特性分析方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
HAO, YAOYAO: "Motor Cortical Representation and Decoding of Monkey Reach and Grasp Movement for Brain Machine Interfaces", CHINA DOCTORAL DISSERTATIONS FULL-TEXT DATABASE (BASIC SCIENCES, 15 August 2013 (2013-08-15), pages 39, ISSN: 1674-022X *

Also Published As

Publication number Publication date
US11238991B2 (en) 2022-02-01
US20190378622A1 (en) 2019-12-12
CN106529186A (zh) 2017-03-22
CN106529186B (zh) 2018-10-02

Similar Documents

Publication Publication Date Title
CN105279365B (zh) 用于学习异常检测的样本的方法
CN112262440B (zh) 一种通过影像组学特征判断癌症治疗反应的方法及系统
WO2018094718A1 (zh) 一种脑区脉冲神经信号的预测方法
WO2019174142A1 (zh) 一种多模式的退化过程建模及剩余寿命预测方法
KR101860061B1 (ko) 심층 신경망 기반 질병 정보 예측 시스템 및 방법
WO2018076571A1 (zh) Lte网络中的异常值检测方法及系统
Mahmud et al. A poisson process model for activity forecasting
EP3629201B1 (en) Engineering automation using quantification and visualization of global sensitivities in spatial temporal domain
CN114358172B (zh) 核反应堆故障分类方法、装置、计算机设备以及存储介质
CN105975772B (zh) 基于概率假设密度滤波的多目标检测前跟踪方法
CN110782181A (zh) 一种低压台区线损率的计算方法及可读存储介质
CN119066913A (zh) 一种适用于多尺度仿真的高性能计算方法
CN116728783A (zh) 一种基于3d打印机的仿真方法及系统
CN114387545A (zh) 基于前馈网络的角膜生物力学特性的智能检测方法
CN117943892A (zh) 基于时序数据分析的刀具监测数据预测优化方法
CN117610725A (zh) 基于多传感数据融合的船舶典型设备寿命预测方法
CN119691410A (zh) 一种基于风险感知的自适应时序预测方法
CN110517774A (zh) 一种预测体温异常的方法
Sun et al. Combined feature selection and cancer prognosis using support vector machine regression
CN114488297B (zh) 一种断层识别方法及装置
CN114491845A (zh) 融合历史轨迹的船用动力轴承剩余寿命预测方法及系统
Jansson et al. Non-parametric analysis of eye-tracking data by anomaly detection
JP4299508B2 (ja) 製造プロセスにおける操業と品質の関連分析装置、関連分析方法及びコンピュータ読み取り可能な記憶媒体
CN104462784A (zh) 一种基于动态分辨熵的传感器优化管理方法
CN109065168B (zh) 一种基于时空聚类统计进行疾病风险评估的方法

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 16922069

Country of ref document: EP

Kind code of ref document: A1

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 16922069

Country of ref document: EP

Kind code of ref document: A1