WO2023040521A1 - 一种自适应闭环深部脑刺激方法、装置及电子设备 - Google Patents
一种自适应闭环深部脑刺激方法、装置及电子设备 Download PDFInfo
- Publication number
- WO2023040521A1 WO2023040521A1 PCT/CN2022/111364 CN2022111364W WO2023040521A1 WO 2023040521 A1 WO2023040521 A1 WO 2023040521A1 CN 2022111364 W CN2022111364 W CN 2022111364W WO 2023040521 A1 WO2023040521 A1 WO 2023040521A1
- Authority
- WO
- WIPO (PCT)
- Prior art keywords
- stimulation
- response
- target
- time
- brain stimulation
- Prior art date
Links
- 230000000638 stimulation Effects 0.000 title claims abstract description 226
- 210000004556 brain Anatomy 0.000 title claims abstract description 219
- 238000000034 method Methods 0.000 title claims abstract description 86
- 230000003044 adaptive effect Effects 0.000 title claims abstract description 32
- 239000002245 particle Substances 0.000 claims abstract description 123
- 230000001537 neural effect Effects 0.000 claims abstract description 26
- 238000004422 calculation algorithm Methods 0.000 claims abstract description 23
- 238000005457 optimization Methods 0.000 claims abstract description 20
- 238000012937 correction Methods 0.000 claims abstract description 8
- 230000004044 response Effects 0.000 claims description 156
- 238000012549 training Methods 0.000 claims description 55
- 208000002193 Pain Diseases 0.000 claims description 51
- 230000000694 effects Effects 0.000 claims description 26
- 238000012360 testing method Methods 0.000 claims description 26
- 230000002123 temporal effect Effects 0.000 claims description 19
- 238000007781 pre-processing Methods 0.000 claims description 17
- 238000004364 calculation method Methods 0.000 claims description 13
- 230000003595 spectral effect Effects 0.000 claims description 13
- 238000010276 construction Methods 0.000 claims description 12
- 238000011156 evaluation Methods 0.000 claims description 11
- 230000008569 process Effects 0.000 claims description 11
- 230000001419 dependent effect Effects 0.000 claims description 8
- 238000000611 regression analysis Methods 0.000 claims description 7
- 238000012545 processing Methods 0.000 claims description 6
- 238000005070 sampling Methods 0.000 claims description 6
- 238000001914 filtration Methods 0.000 claims description 5
- 239000000284 extract Substances 0.000 claims description 4
- 238000012544 monitoring process Methods 0.000 claims description 4
- 238000012847 principal component analysis method Methods 0.000 claims description 3
- 230000010354 integration Effects 0.000 claims 1
- 230000007774 longterm Effects 0.000 abstract description 6
- 230000000875 corresponding effect Effects 0.000 description 67
- 230000006870 function Effects 0.000 description 8
- 238000010586 diagram Methods 0.000 description 7
- 238000001356 surgical procedure Methods 0.000 description 7
- 230000008859 change Effects 0.000 description 6
- 238000010219 correlation analysis Methods 0.000 description 6
- 208000004296 neuralgia Diseases 0.000 description 5
- 208000021722 neuropathic pain Diseases 0.000 description 5
- 230000001133 acceleration Effects 0.000 description 3
- 230000009471 action Effects 0.000 description 3
- 230000002596 correlated effect Effects 0.000 description 3
- 238000012804 iterative process Methods 0.000 description 3
- 230000004007 neuromodulation Effects 0.000 description 3
- 230000033764 rhythmic process Effects 0.000 description 3
- 238000001228 spectrum Methods 0.000 description 3
- 238000004458 analytical method Methods 0.000 description 2
- 230000009286 beneficial effect Effects 0.000 description 2
- 239000000090 biomarker Substances 0.000 description 2
- 230000007177 brain activity Effects 0.000 description 2
- 238000001514 detection method Methods 0.000 description 2
- 229940079593 drug Drugs 0.000 description 2
- 239000003814 drug Substances 0.000 description 2
- 238000000605 extraction Methods 0.000 description 2
- 230000004048 modification Effects 0.000 description 2
- 238000012986 modification Methods 0.000 description 2
- 210000004761 scalp Anatomy 0.000 description 2
- 238000010183 spectrum analysis Methods 0.000 description 2
- 230000004936 stimulating effect Effects 0.000 description 2
- NAWXUBYGYWOOIX-SFHVURJKSA-N (2s)-2-[[4-[2-(2,4-diaminoquinazolin-6-yl)ethyl]benzoyl]amino]-4-methylidenepentanedioic acid Chemical compound C1=CC2=NC(N)=NC(N)=C2C=C1CCC1=CC=C(C(=O)N[C@@H](CC(=C)C(O)=O)C(O)=O)C=C1 NAWXUBYGYWOOIX-SFHVURJKSA-N 0.000 description 1
- 208000012902 Nervous system disease Diseases 0.000 description 1
- 208000025966 Neurological disease Diseases 0.000 description 1
- 208000004550 Postoperative Pain Diseases 0.000 description 1
- 230000002159 abnormal effect Effects 0.000 description 1
- 230000002411 adverse Effects 0.000 description 1
- 230000004888 barrier function Effects 0.000 description 1
- 230000008901 benefit Effects 0.000 description 1
- 210000003710 cerebral cortex Anatomy 0.000 description 1
- 208000010877 cognitive disease Diseases 0.000 description 1
- 238000013135 deep learning Methods 0.000 description 1
- 238000013136 deep learning model Methods 0.000 description 1
- 238000010828 elution Methods 0.000 description 1
- 238000005516 engineering process Methods 0.000 description 1
- 230000019637 foraging behavior Effects 0.000 description 1
- 230000006872 improvement Effects 0.000 description 1
- 238000012417 linear regression Methods 0.000 description 1
- 230000000926 neurological effect Effects 0.000 description 1
- 210000002569 neuron Anatomy 0.000 description 1
- 230000001575 pathological effect Effects 0.000 description 1
- 230000035790 physiological processes and functions Effects 0.000 description 1
- 230000002980 postoperative effect Effects 0.000 description 1
- 238000000513 principal component analysis Methods 0.000 description 1
- 238000003672 processing method Methods 0.000 description 1
- 208000020016 psychiatric disease Diseases 0.000 description 1
- 238000012216 screening Methods 0.000 description 1
- 238000010845 search algorithm Methods 0.000 description 1
- 210000003625 skull Anatomy 0.000 description 1
- 238000012421 spiking Methods 0.000 description 1
- 230000003068 static effect Effects 0.000 description 1
- 238000007619 statistical method Methods 0.000 description 1
- 230000001960 triggered effect Effects 0.000 description 1
- 238000010200 validation analysis Methods 0.000 description 1
Images
Classifications
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
- A61B5/00—Measuring for diagnostic purposes; Identification of persons
- A61B5/24—Detecting, measuring or recording bioelectric or biomagnetic signals of the body or parts thereof
- A61B5/316—Modalities, i.e. specific diagnostic methods
- A61B5/388—Nerve conduction study, e.g. detecting action potential of peripheral nerves
-
- A—HUMAN NECESSITIES
- A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
- A61N—ELECTROTHERAPY; MAGNETOTHERAPY; RADIATION THERAPY; ULTRASOUND THERAPY
- A61N1/00—Electrotherapy; Circuits therefor
- A61N1/18—Applying electric currents by contact electrodes
- A61N1/32—Applying electric currents by contact electrodes alternating or intermittent currents
- A61N1/36—Applying electric currents by contact electrodes alternating or intermittent currents for stimulation
-
- 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
Definitions
- the invention relates to the field of medical electronic systems, in particular to an adaptive closed-loop deep brain stimulation method, device and electronic equipment.
- DBS Deep brain stimulation
- the continuous open-loop stimulation mode is used clinically, and the doctor adjusts the stimulation parameters according to the patient's condition and fixes them until the next visit.
- Stimulation parameters include stimulation amplitude, frequency, and pulse width.
- the stimulation parameters cannot be properly adjusted according to the patient's instantaneous or long-term state changes, and long-term continuous stimulation may also bring about many side effects such as language barriers and cognitive dysfunction.
- the closed-loop DBS system controls the neural activity of the brain by applying electrical stimulation signals that can be adjusted according to the control target. How to form an adaptive closed-loop DBS by detecting pathological neural activity has become an important problem to be solved to improve the clinical treatment effect of DBS.
- the purpose of the present invention is to provide an adaptive closed-loop deep brain stimulation method, device and electronic equipment, which can effectively solve the problem of personalized neural regulation under multi-state and long-term conditions.
- an adaptive closed-loop deep brain stimulation method comprising:
- the parameters of the target proportional-derivative-integral controller are corrected online while performing deep brain stimulation.
- the parameter search through the particle swarm optimization algorithm to obtain the target proportional-derivative-integral controller parameters to determine the target proportional-derivative-integral controller includes:
- the position coordinates of the particles corresponding to the global optimal fitness after the current iteration are used as the parameters of the target proportional-derivative-integral controller to determine the target proportional-derivative-integral controller.
- the judging whether the particle swarm after the current iteration meets the condition for terminating the iteration includes:
- the iterative calculation of the current fitness of any particle in the particle swarm within any first window length and updating the global optimal fitness includes:
- the current fitness of any particle in the current first window length is less than its previous fitness, the current fitness of the corresponding particle is taken as the individual optimal fitness of the corresponding particle;
- the individual optimal fitness of any particle is smaller than the individual optimal fitness of other particles in the particle swarm, the individual optimal fitness of the corresponding particle is taken as the global optimal fitness.
- the iterative calculation of the current fitness of any particle in the particle swarm within any first window length includes:
- the current fitness is obtained based on the current neural activity signal and a preset target signal.
- the deep brain stimulation using the stimulation parameters obtained by the target proportional-derivative-integral controller includes:
- Stimulation pulses are formed based on the target stimulation parameters to perform deep brain stimulation.
- the online correction of the target proportional-derivative-integral controller parameters while performing deep brain stimulation includes:
- the monitoring of neural activity signals during deep brain stimulation and judging whether to adjust the target proportional-derivative-integral controller parameters include:
- the method further includes: acquiring a target brain stimulation response corresponding to adaptive closed-loop deep brain stimulation, including:
- a target brain stimulus response is generated through a pre-built brain stimulus response model.
- the method also includes:
- the brain stimulation response model is constructed in advance, and the construction method includes:
- the training sample set including at least one set of temporal stimulus inputs to the brain and a temporal real response corresponding to each set of temporal stimulus inputs;
- the brain stimulation response model is obtained based on generative confrontation network training.
- the confrontation network includes a generation network and a confrontation network
- Said taking said timing stimulus input and said timing real response as input, and obtaining said brain stimulation response model based on generative confrontation network training include:
- timing stimulus input as the input of the generation network to obtain a corresponding timing generation response
- the training is stopped and the model corresponding to the generation network is used as the brain stimulation response model.
- the training is stopped and the model corresponding to the generation network is used as the brain stimulation response model, including:
- the said brain stimulation response model is obtained based on the training of generative confrontation network by taking said time series stimulus input and said time series real response as input, further comprising:
- the backpropagation algorithm is used to update the weights and biases of the generation network and the confrontation network.
- the method further includes: performing model evaluation on the brain stimulation response model based on a pre-acquired test sample set, including:
- the time-sequence stimulation input includes a time-sequence stimulation amplitude and a time-sequence stimulation frequency
- the time-sequence real response includes a time-series real local field potential corresponding to the time-sequence stimulation amplitude and time-sequence stimulation frequency collected signal, the temporally generated response comprising a temporally generated local field potential signal generated by the generating network based on the stimulus amplitude and the stimulus frequency.
- the method before obtaining the brain stimulation-response model based on Generative Adversarial Network training using the time-sequence stimulus input and the time-sequence real response as input, the method further includes analyzing the collected The time series real response is preprocessed, including:
- the corresponding power time series is calculated to obtain the preprocessed time-series real local field potential signal.
- the method also includes: constructing a pain state prediction model based on brain electrical signals, including:
- the characteristics of brain electrical activity are screened in the time domain and wavelet domain; in the frequency domain, the principal component analysis method PCA is used to obtain the key components that characterize each feature group according to the contribution rate, and then the brain electrical activity is screened. activity characteristics;
- the features extracted in the time domain include the average value, standard deviation and information entropy of the signal amplitude; the signal is standardized before the feature is extracted, and the specific method is to divide the signal value of each sampling point by the maximum amplitude value; the features extracted by the wavelet domain are the percentage of the synchronization state existence time in the delta, theta, alpha, low-beta, high-beta, low-gamma and high-gamma frequency bands and the 21 The percentage of the occurrence time of the four states 00, 01, 10, and 11 composed of the binary codes of the synchronization levels of each frequency segment in the combination state to the total time.
- step 1) the preprocessing of the brain electrical signal includes:
- the feature extracted in the frequency domain is the power value of the power spectral density integrated on different frequency segments after the Fourier transform and the ratio of power between different frequency bands; before extracting the feature
- the signal is standardized, and the specific method is to divide the power spectral density value at each frequency point by the integral of the power spectral density in the frequency range of 2-90Hz.
- step 3 in the time domain and wavelet domain, select the feature with the significance of the pain state less than 0.05 or 0.01; Component selection features.
- step 5 every time the pain state prediction model completes a prediction, current data is incorporated into it to modify model parameters.
- an adaptive closed-loop deep brain stimulation device comprising:
- a parameter search module configured to perform parameter search through a particle swarm optimization algorithm to obtain target proportional-derivative-integral controller parameters to determine the target proportional-derivative-integral controller;
- a stimulation module for performing deep brain stimulation using stimulation parameters obtained by the target proportional-derivative-integral controller
- the correction module is used for performing online correction on the target proportional-derivative-integral controller parameters while performing deep brain stimulation.
- the device also includes:
- the first acquisition module is used to acquire target timing stimulation input
- a generating module configured to generate a target brain stimulus response through a pre-built brain stimulus response model based on the target timing stimulus input.
- an electronic device including:
- a memory associated with the one or more processors the memory is used to store program instructions, and when the program instructions are read and executed by the one or more processors, the execution of any one of the first aspect described method.
- the application provides an adaptive closed-loop deep brain stimulation method, device and electronic equipment, wherein the method includes: performing parameter search through particle swarm optimization algorithm to obtain target proportional-differential-integral controller parameters, and adopting the target proportional-differential-integral control
- the stimulation parameters corresponding to the device parameters are used for deep brain stimulation, and the target proportional-derivative-integral controller parameters are corrected online while deep brain stimulation is performed.
- This method can automatically calculate the controller gain for different patients and change with the patient state. Automatic calibration of controller gains for personalized neuromodulation under multi-state, long-term conditions;
- the application also includes obtaining the target brain stimulation response corresponding to the adaptive closed-loop deep brain stimulation, including obtaining the target timing stimulation input; based on the target timing stimulation input, generating the target brain stimulation response through a pre-built brain stimulation response model; the method According to the time-varying, nonlinear and uncertain characteristics of the brain in the process of responding to stimulation, it can truly simulate the response of the brain after receiving deep brain stimulation, so as to improve the accuracy of brain stimulation parameters;
- the application also includes the construction of a pain state prediction model based on brain electrical signals, which can comprehensively characterize and quantify brain electrical activity from a multi-dimensional perspective, combine subjective evaluation and objective detection methods, and combine the detected multiple Biomarkers are fused to establish a quantitative prediction model of the patient's state or degree of change, which can be used for accurate judgment of pain or pain status;
- Fig. 1 is the flowchart of adaptive closed-loop deep brain stimulation method in the present embodiment
- Fig. 2 is another flowchart of the self-adaptive closed-loop deep brain stimulation method in this embodiment
- FIG. 3 is a schematic diagram of parameter calculation and adaptive stimulation timing in this embodiment
- FIG. 4 is a schematic diagram of the closed-loop deep brain stimulation principle based on a PID controller in this embodiment
- Fig. 5 is the flow chart of the method for constructing the brain stimulation response model in the present embodiment
- Fig. 6 is a schematic diagram of a method for constructing a brain stimulation response model in this embodiment
- Fig. 7 is a flow chart of the brain stimulation response method in this embodiment.
- Fig. 8 is a schematic flowchart of an example of constructing a prediction model based on brain electrical signal pain state (pain relief degree of pain patients).
- FIG. 9 is an example diagram of the calculation results of the power of the frequency point and the ratio of the power between frequency bands of the electrical signal (ie, LFP) recorded in the brain of a patient with neuropathic pain using power spectrum analysis.
- Fig. 10 is a diagram of the correlation coefficient result obtained from the correlation analysis between the frequency domain features and the degree of pain relief of the patient, and shows the identified features related to the degree of pain relief of the patient.
- Figure 11 is an example graph of performance comparison and validation results of prediction models obtained in pain patient data.
- first and second are used for description purposes only, and cannot be interpreted as indicating or implying relative importance or implicitly indicating the quantity of indicated technical features. Thus, a feature defined as “first” and “second” may explicitly or implicitly include one or more of these features. In the description of the present invention, unless otherwise specified, "plurality" means two or more.
- the brain has nonlinear, time-varying, and non-stationary characteristics. It is difficult to generalize how to drive individual brain network dynamics by using the traditional spiking neuron model and LSSM model, or it does not conform to the nonlinear characteristics of the brain's response to stimuli. In view of the above-mentioned current situation of brain stimulation, it is necessary to find a brain stimulation response method corresponding to the above-mentioned characteristics of the brain, so as to truly simulate the response of the brain when it is stimulated by the brain.
- the brain stimulation response model construction method, response method, device and electronic equipment of this embodiment will be further described in detail below with reference to FIGS. 1 to 11 .
- this embodiment provides an adaptive closed-loop deep brain stimulation method, including the following steps:
- PSO is a random search algorithm based on group cooperation developed by simulating the foraging behavior of birds, which can be used to solve optimization problems.
- Each bird is abstracted as a "particle" in the parameter space, and the position of the food is abstracted as a feasible solution that meets the requirements.
- All particles have a fitness value determined by the function being optimized, and the velocity of each particle determines the direction and distance of their parameter search, which is updated after each generation of calculation.
- the particles search in the solution space guided by the global optimal solution and the individual historical optimal solution.
- the proportional-differential-integral (Proportional-Integral-Differential, PID) controller includes a total of three control links: a proportional link, an integral link, and a differential link.
- the integral link produces a control effect based on the historical error, which is mainly used to eliminate the static error and improve the accuracy of the system.
- the strength of the integral action depends on the integral time constant. The larger the gain of the integral link, the weaker the integral action, and vice versa.
- the differential link can reflect the change trend (change rate) of the deviation signal, and can introduce an effective early correction signal into the system before the deviation signal value becomes too large, thereby speeding up the action speed of the system and reducing the adjustment time.
- the PID controller has the characteristics of simple algorithm, good robustness and high reliability.
- an incremental PID controller is specifically used in this embodiment, which is a modification of the classic PID.
- step S100 includes:
- the PSO parameters shown in Table 1 are initialized.
- the particle position x and particle velocity v are random values;
- the particle position dimension d corresponds to the number of controller gains to be determined;
- the number of particles N, the maximum number of iterations G, the first window length t 1 of each iteration, and the inertia weight w , acceleration constant c 1 /c 2 , and initial value range b are artificially set values.
- the value range of the initial position range b will affect the time-consuming of the PSO search process. If there is no limit, since the initial position is random, it may be far away from the optimal solution that meets the requirements, and the search process will take a long time.
- step S12 Iteratively calculate the current fitness of any particle in the particle swarm within any first window length and update the global optimal fitness.
- step S12 includes:
- the PID controller is used to perform deep brain stimulation with a certain gain and obtain neural activity signals while stimulating. After the signal is processed, it obtains a fitness (steady-state error) compared with the preset target signal, which is used to adjust the gain of the PID controller. Adjustments are made to achieve real-time closed-loop control.
- S121 specifically includes:
- the current position coordinate is the particle position vector.
- the stimulation parameters were calculated using the following formulas (1) and (2):
- u(k) is the stimulus parameter.
- u(k) is at least one of amplitude, frequency, and pulse width.
- the frequency and pulse width are preferably set values, such as 130 Hz and 60 ⁇ s, so u(k) is the amplitude.
- the brain stimulation amplitude is obtained through the above formulas (1) (2), deep brain stimulation is performed, and the stimulation duration is the first window length t 1 .
- each of the N first window lengths t1 needs to be sequentially performed with the current position coordinates of each random particle as the PID gain for deep brain stimulation.
- the current position coordinates of the random particle 1 in the first window length t1 as the PID gain to perform deep brain stimulation and calculate the current position coordinates of the particle 1 in the first window length t1 adaptability.
- the current coordinates of the particle 2 are used as the PID gain to perform deep brain stimulation, as a.
- an elution period t2 as shown in FIG. 4 is set.
- a washout period t2 is set after each first window length t1 , and the brain is not stimulated in any way during the washout period t2 , so as to ensure that each step in the particle swarm optimization process Brain states consistent with deep brain stimulation.
- the corresponding fitness is calculated through a fitness function.
- the control target as an example to suppress the energy of a certain frequency band of neural activity
- the fitness function is to minimize the energy of this frequency band and the preset target
- t s represents the start time of each time window
- t e represents the end time of each time window
- e(t) represents the error at each feedback moment in the current window length.
- control target can be expanded from a single target to multiple targets, and when extended to multi-target control, the dimension of the particle position vector changes accordingly.
- the fitness function can be modified as formula (4) :
- w 1 and w 2 represent the inertia weight of error and control quantity to parameter selection, respectively.
- this method may be limited when the number of control targets is too large.
- the current fitness of any particle in the current first window length is less than its previous fitness
- the current fitness of the corresponding particle is taken as the individual optimal fitness of its individual optimal position pbest, that is, carry out The individual optimal position pbest of the single particle itself and the update of the individual optimal fitness.
- the particle position of the particle is the global best position gbest.
- the particle velocity and particle position can be updated, as shown in formulas (5) and (6):
- v i ⁇ v i +c 1 rand(0,1)(pbest i,d -xi )+c 2 rand(0,1)(gbest d -xi ) (5)
- pbest i,d is the individual optimal position of the i-th particle at the d-th iteration
- gbest d is the global optimal position at the d-th iteration
- w is the inertia weight
- c 1 and c 2 are the acceleration constants respectively.
- S130 in the above update process, judge whether the particle swarm after the current iteration meets the termination iteration condition; if so, use the position coordinates of the particles corresponding to the global optimal fitness after the current iteration as the target PID controller parameter, thereby Determine the target proportional-derivative-integral controller.
- the conditions for terminating the iteration including but not limited to the judgment result of any one of the number of iterations, the average fitness, or the global optimal fitness.
- step S130 includes:
- the preferred frequency and pulse width are preferably set values, so the target stimulation parameter u(k) refers to the amplitude.
- S220 forming stimulation pulses based on the target stimulation parameters to perform deep brain stimulation.
- the adaptive closed-loop deep brain stimulation method also includes:
- Step S30 Perform online calibration of the target PID controller parameters while performing deep brain stimulation.
- Step S30 specifically includes:
- S310 Monitor the neural activity signal during the deep brain stimulation and judge whether it is necessary to adjust the target PID controller parameters; including:
- the steady-state error is the aforementioned fitness. Since the physiological state of the stimulated object is changing in real time, especially when the condition of the stimulated object changes, medication, exercise and other states change, the fitness under the same stimulation parameters will increase significantly, and continuous stimulation will cause serious damage to the stimulated object. adverse effects.
- step S200 If not, continue to execute step S200.
- the adaptive closed-loop deep brain stimulation method provided in this embodiment can automatically calculate the controller gain for different patients and automatically calibrate the controller gain as the patient's state changes, so as to realize personalization under multi-state and long-term conditions neuromodulation.
- the adaptive closed-loop deep brain stimulation method in this embodiment further includes S400.
- step S400 includes the following steps:
- the stimulation input is recorded, while the real response output of the brain is collected.
- the stimulus input includes stimulus amplitude U and stimulus frequency f
- the real response output is a local field potential (local field power, LFP) signal y.
- the data in the training sample set are given time sequence, that is, the stimulus input is time sequence stimulus input x k
- the real output response is also time sequence real response y k .
- the temporal stimulus input x k includes the temporal stimulus amplitude and temporal stimulus frequency
- the temporal real response y k includes the collected temporal real partial field potential signal corresponding to the temporal stimulus amplitude and temporal stimulus frequency
- step S420 it also includes: S40, preprocessing the collected time series real response y k , specifically including:
- the broadband original signal recorded during the stimulation process contains stimulation artifacts.
- the template method is a commonly used stimulation artifact removal method in the field of brain stimulation, which will not be further described in this embodiment.
- an anti-aliasing filter with a cutoff frequency of 100 Hz is used to down-sample the time-series real local field potential signal to 200 Hz.
- this embodiment adopts the use of passband cutoff frequency 1Hz, the equiripple finite impulse response (FIR) filter of stopband cutoff frequency 0.5Hz to remove drift, uses stopband cutoff frequency 59Hz and 61Hz and passband cutoff frequency 58Hz and 62Hz band-stop equiripple FIR filters to remove line noise at 60Hz, and band-stop equiripple FIR filters with stop-band cut-off frequencies of 49Hz and 51Hz and pass-band cut-off frequencies of 48Hz and 52Hz to eliminate any possible residual at the stimulus frequency stimulus artifacts.
- FIR finite impulse response
- each LFP channel is divided into multiple time windows in sequence according to the preset time window length T w and the average power of the LFP in each time window is calculated, thus obtaining the power time series y k of the LFP, that is, the sequence True local field potential signal.
- the timing stimulus input x k may be a preset empirical value.
- the selectable value of the temporal stimulation amplitude is 0V (no stimulation), 1.5V or 3V
- step S2 After determining the training sample set and performing corresponding data preprocessing, the model training in step S2 is performed, and the specific training is as follows.
- the generative adversarial network is a deep learning model, which includes two network models, the generative network and the discriminant network, which are set in sequence.
- the task of generative network is to generate examples that look natural and real, similar to the original data.
- the task of discriminative networks is to judge whether a given instance looks natural or artificial.
- Generative adversarial networks are trained to simulate real brain responses to brain stimuli by competing with a non-linear generator that generates responses and a discriminator that discriminates whether responses are real or fake.
- the input of the discriminative network is the sequence generation response produced by the generative network and the temporal real response of the brain to stimuli y k .
- step S420 specifically includes:
- step S423 includes:
- timing generation response is the same as the temporal real response y k , then stop the training and use the model corresponding to the generation network as the brain stimulus response model; or,
- the judgment result is timing generation response Different from the time series real response y k , continue training until the judgment result meets the preset threshold, stop training and use the model corresponding to the generation network as the brain stimulation response model. In this step, it is used to generate the response when the timing When it is different from the real response y k of the time series, evaluate the closeness of the two.
- the preset threshold is preferably 0.5, that is, when the judgment result is closer to 0.5, the timing generated response The closer it is to the real response y k of time series, the higher the authenticity.
- a backpropagation algorithm is used to update the weights and biases of the generation network and the confrontation network.
- G(x) is the generation network
- D(x) is the confrontation network
- the generated network update gradient is shown in the following formula (9):
- the method further includes: S430, performing model evaluation on the brain stimulation response model based on the pre-acquired test sample set, specifically including:
- the test sample set includes time series stimulus input x k and the corresponding time series true response y k , and time series stimulus input x k includes time series stimulus amplitude and time series stimulus frequency.
- CC Pearson's correlation coefficient
- CC is used to measure the degree of correlation between two variables, and its value is between -1 and 1.
- CC is based on the timing of GAN to generate responses The amount of the linear correlation degree between the time series real response y k , the larger the value, the higher the correlation degree, and the higher the accuracy of the brain stimulation response model.
- CC expression in the present embodiment is as shown in formula (10):
- Cov() and Var() represent the covariance and variance of the time series, respectively.
- step S400 also includes S440 obtaining the target brain stimulation response corresponding to the adaptive closed-loop deep brain stimulation.
- the brain stimulation response method includes:
- the above-mentioned target timing stimulus input includes target stimulus amplitude U and target stimulus frequency f.
- the target brain stimulus response corresponding to the target timing stimulus input obtained through the brain stimulus response model can effectively simulate the real LFP signal.
- the brain stimulation response model construction method provided in this embodiment is based on the generation of adversarial network modeling to obtain the obtained brain stimulation response model. Based on the powerful learning ability of deep learning, the brain stimulation response model can target the The characteristics of time-varying, non-linear and uncertain, etc., truly simulate the stimulus response of the brain after being stimulated.
- the adaptive closed-loop deep brain stimulation method provided in this embodiment also includes building a pain state prediction model based on brain electrical signals, as shown in Figure 8, including:
- the signal can be the scalp electroencephalogram (EEG) recorded by the scalp electrode, the electrocortical electroencephalogram (ECoG) signal placed under the skull and in contact with the cerebral cortex, or the signal implanted into the brain.
- EEG scalp electroencephalogram
- EoG electrocortical electroencephalogram
- LFPs Field potential signals of deep nuclei associated with pain.
- the recording time is not less than 180s. It is also necessary to record the results of the assessment of the pain patient's status by the clinical staff.
- a common assessment scale is the Pain Analogue Scale (VAS).
- the features extracted in the time domain include the mean value, standard deviation and information entropy of the signal amplitude. Signals are normalized before feature extraction. The specific method is to divide the signal value of each sampling point by the maximum value of the amplitude.
- the frequency domain feature is the power value of the power spectral density integrated in different frequency bands after Fourier transform and the ratio of power between different frequency bands.
- the integral of power spectral density in a certain frequency interval represents the activity level of the signal in this frequency band, and the ratio of activity levels between different frequency bands is also used as a feature of brain activity.
- the power spectrum of each patient is standardized, and the processing method is to divide the power spectral density value at each frequency point by the integral of the power spectral density in the frequency range of 2-90Hz; different
- the ratio of power between frequency bands refers to the ratio between the amplitude values of signals in different frequency bands.
- the method of traversing the frequency combination is adopted. Considering that the commonly used rhythm frequency band is generally 4Hz or its multiples , so the analysis uses a frequency band of 4 Hz with a step size of 0.5 Hz to traverse the frequency combinations.
- the wavelet domain feature adopts the method designed in the applicant's patent 201610487800.X to discriminate the synchronization status of signals in 7 frequency bands. These frequency bands are delta (3-6Hz), theta (6-9Hz), alpha (9-12Hz), low-beta (12-24Hz), high-beta (24-36Hz), low-gamma (36-60Hz ), high-gamma (60-90Hz). And the 7 frequency band signals are combined in pairs to obtain 20 combination states, and each combination has four states, which are 00, 01, 10, and 11 respectively. Calculate the percentage of the total time of the single frequency over-synchronization state and the combined state of the four states of the total time respectively.
- the synchronization time of only the highbeta of a single rhythm was significantly correlated with the degree of pain relief at 12 months after surgery.
- Multiple states in the combined state consisting of both rhythms were associated with greater pain relief at 12 months postoperatively.
- Table 1 shows the results of the correlation analysis between the time-domain features and the degree of pain and pain relief 12 months after surgery in patients with neuropathic pain.
- Table 2 shows the results of the correlation analysis between the percentage of time in the state of excessive synchronization in the 7 frequency bands of wavelet domain characteristics and the degree of pain and pain relief in patients with neuropathic pain at 12 months after surgery.
- Table 3 is a partial result showing the significant correlation between the percentage of the synchronization state existence time of the combination of 7 frequency bands of the wavelet domain feature and the degree of pain relief 12 months after surgery in patients with neuropathic pain.
- FIG. 9 is an example diagram of the calculation results of the power of the frequency point and the power between frequency bands of the electrical signal (ie, LFP) recorded in the brain of a patient with neuropathic pain using power spectrum analysis.
- Fig. 10 is a diagram of the correlation coefficient result obtained from the correlation analysis between the frequency domain features and the degree of pain relief of the patient, and shows the identified features related to the degree of pain relief of the patient.
- the power spectrum characteristics of the frequency of the thick solid line in the left figure of Figure 10 are significantly correlated with the degree of pain relief at 12 months after surgery (P ⁇ 0.01), while the power spectrum between the frequency bands in the boxed part in the right figure of Figure 10
- the ratio characteristic of density was significantly correlated with the degree of pain relief 12 months after operation (P ⁇ 0.01).
- the method of using principal component analysis in modeling only selects the first principal component in Fig. 9 and Fig. 10 .
- the characteristics selected from the three dimensions are used as independent variables, and the pain degree or degree of relief at 12 months after operation is used as the dependent variable, and the prediction models are respectively constructed through regression analysis.
- the result of the prediction model was used as the independent variable, and the degree of pain or pain relief at 12 months after operation was used as the dependent variable.
- a pain state prediction model integrating three-dimensional features was established through multiple linear regression analysis.
- Figure 11 shows the specific prediction results. The results show that the predictive effect of integrating the features of time domain, frequency domain and wavelet domain is better than that of only using single dimension features.
- the integrated prediction of pain severity at 12 months after surgery can reach 75%, and the integrated prediction accuracy of pain relief at 12 months after surgery can reach 83%.
- this embodiment comprehensively characterizes and quantifies brain electrical activity from a multi-dimensional perspective, combines subjective evaluation and objective detection means, and separates and fuses multiple biomarkers from a single brain electrical activity to achieve state Accurate judgment of patient status and quantitative prediction model of change degree.
- the pain state can be effectively predicted by the model of the present invention.
- this embodiment further provides an adaptive closed-loop deep brain stimulation device, including:
- the parameter search module is used for performing parameter search through the particle swarm optimization algorithm to obtain target proportional-differential-integral controller parameters to determine the target proportional-differential-integral controller. Further, the parameter search module includes:
- the initialization unit is used to initialize the parameters of the particle swarm optimization algorithm.
- the update unit is used to iteratively calculate the current fitness of any particle in the particle swarm within any first window length and update the global optimal fitness.
- the update unit includes:
- the first calculation subunit is used to iteratively calculate the current fitness of any particle in the particle swarm within any first window length.
- the first calculation subunit is used for:
- the current fitness is obtained based on the current neural activity signal and a preset target signal.
- the first update subunit is used to take the current fitness of the corresponding particle as the individual optimal fitness of the corresponding particle when the current fitness of any particle in the current first window length is less than any previous fitness Spend.
- the second update subunit is used to take the individual optimal fitness of the corresponding particle as the global optimal when the individual optimal fitness of any particle is smaller than the individual optimal fitness of the rest of the particles in the particle swarm adaptability.
- the first judging unit is used to judge whether the particle swarm after the current iteration meets the condition for terminating the iteration. Specifically, the first judging unit is used for:
- the first update unit is used to determine the target proportional-differential-integral controller parameters by using the position coordinates of the particles corresponding to the global optimal fitness after the current iteration as the target proportional-differential-integral controller parameter when the judgment result is yes. Integral controller.
- the stimulus module includes:
- the stimulation unit is used to form stimulation pulses based on the target stimulation parameters to perform deep brain stimulation.
- the device also includes a correction module, which is used for online correction of the target proportional-derivative-integral controller parameters while performing deep brain stimulation.
- the calibration module includes:
- the second judging unit is used to monitor the neural activity signal during the deep brain stimulation and judge whether it is necessary to adjust the target proportional-derivative-integral controller parameters;
- the second judging unit includes:
- the second calculation subunit is used to obtain the corresponding steady-state error based on the neural activity signal and the preset target signal during deep brain stimulation in any second window length;
- a judging subunit configured to judge that the target proportional-derivative-integral controller parameter needs to be adjusted when the number of consecutive second window lengths in which the steady-state error exceeds a preset steady-state error threshold reaches a preset window number.
- the adaptive closed-loop deep brain stimulation device also includes:
- a first acquisition module configured to acquire a training sample set, the training sample set including at least one set of time-sequence stimulation inputs to the brain and a time-series true response corresponding to each set of time-sequence stimulation inputs;
- the first training module is used to obtain the brain stimulation response model based on the training of the generative confrontation network by taking the time series stimulus input and the time series real response as input.
- the generative confrontation network includes a generative network and a confrontation network.
- the first training module includes:
- a generation network model unit is used to use the timing stimulus input as the input of the generation network to obtain a corresponding timing generation response
- the confrontation network model unit is used to use the timing real response corresponding to the timing stimulus input and the timing generated response as the input of the confrontation network to obtain a corresponding judgment result; when the judgment result meets the preset condition , stop training and use the model corresponding to the generation network as the brain stimulus response model.
- the second update unit is used to update the weights and biases of the generation network and the confrontation network by using the backpropagation algorithm during training.
- confrontation network model unit is specifically used for:
- the evaluation module is used to perform model evaluation on the brain stimulation response model based on the pre-acquired test sample set, including:
- test unit configured to input any time-series stimulus input in the test sample set into the brain stimulus response model and obtain corresponding test results
- the processing unit is configured to calculate the Pearson correlation coefficient between the time-series real response corresponding to the time-sequence stimulus input in the test sample set and the corresponding test result, and when the Pearson correlation coefficient meets a preset threshold, the test is passed.
- the time series stimulation input includes time series stimulation amplitude and time series stimulation frequency
- the time series real response includes the collected time series real local field potential signal corresponding to the time series stimulation amplitude and time series stimulation frequency
- the time series generated response includes Local field potential signals are generated based on the timing sequence generated by the generation network based on the stimulation amplitude and stimulation frequency.
- the device for constructing the brain stimulation response model also includes a preprocessing module, which is used for preprocessing the collected time-series real responses.
- the preprocessing modules include:
- a first preprocessing unit configured to remove stimulus artifacts from the collected time-series real local field potential signals
- the second preprocessing unit performs down-sampling on the time-series real local field potential signal from which stimulus artifacts have been removed;
- the third preprocessing unit performs filtering processing on the time-series real local field potential signal after downsampling
- the fourth preprocessing unit calculates a corresponding power time series based on the filtered time-series real local field potential signal to obtain a preprocessed time-series real local field potential signal.
- the unit also includes:
- the second acquisition module is used to acquire target timing stimulation input
- a generating module configured to generate a target brain stimulus response through a brain stimulus response model based on the target timing stimulus input.
- the adaptive closed-loop deep brain stimulation device also includes:
- the second training module including:
- the fifth preprocessing unit is used to preprocess the electrical signals of the brain, and remove signals with poor quality and noises in the signals.
- the preprocessing includes the step of removing 50Hz power frequency interference and baseline drift, and the step of normalizing the noise-removed signal.
- the representation unit is used to extract features from three dimensions of time domain, frequency domain and wavelet domain to characterize brain electrical activity.
- the features extracted in the time domain include the average value, standard deviation, and information entropy of the signal amplitude; the signal is standardized before the feature is extracted, and the specific method is to divide the signal value of each sampling point by the maximum value of the amplitude;
- the features extracted in the wavelet domain are the percentages of the synchronization state existence time in the delta, theta, alpha, low-beta, high-beta, low-gamma and high-gamma frequency bands and 21 combinations obtained by combining these 7 frequency bands in pairs
- the percentage of occurrence time of the four states 00, 01, 10, and 11 composed of the binary codes of the synchronization level of each frequency segment in the state to the total time.
- the feature extracted in the frequency domain is the power value of the power spectral density integrated in different frequency bands after Fourier transform and the ratio of the power between different frequency bands; before the feature is extracted, the signal is standardized.
- the power spectral density value is divided by the integral of the power spectral density in the 2-90Hz frequency band.
- the screening unit is used to screen the characteristics of brain electrical activity in the time domain and wavelet domain according to the correlation with the pain state; in the frequency domain, the principal component analysis method PCA is used to obtain the key components that characterize each feature group according to the contribution rate and then screen Characteristics of brain electrical activity.
- the time domain and wavelet domain select the features whose pain state significance is less than 0.05 or 0.01; in the frequency domain, select the features according to the 1-3 principal components with the largest contribution rate.
- the model construction unit is used to use the features screened from the three dimensions of time domain, frequency domain and wavelet domain as independent variables, and the degree of pain relief as dependent variables to establish state prediction models through regression analysis; predict the state in different dimensions
- the result was used as the independent variable, and the patient's clinical subjective assessment was used as the dependent variable.
- Multiple regression analysis was used to establish an integrated pain state prediction model. Every time the pain state prediction model completes a prediction, the current data is incorporated into it to modify the model parameters.
- the adaptive closed-loop deep brain stimulation device when the adaptive closed-loop deep brain stimulation device provided by the above-mentioned embodiment triggers the adaptive closed-loop deep brain stimulation service, it only uses the division of the above-mentioned functional modules as an example. In practical applications, the above-mentioned Function allocation is accomplished by different functional modules, that is, the internal structure of the system is divided into different functional modules to complete all or part of the functions described above. In addition, the adaptive closed-loop deep brain stimulation device and the adaptive closed-loop deep brain stimulation method provided by the above-mentioned embodiments belong to the same idea, that is, the system is based on this method, and its specific implementation process is detailed in the method embodiment, which is not described here. Let me repeat.
- this embodiment also provides an electronic device, including:
- a memory associated with the one or more processors the memory is used to store program instructions, and when the program instructions are read and executed by the one or more processors, perform the aforementioned adaptive closed-loop deep brain stimulation method.
Landscapes
- Health & Medical Sciences (AREA)
- Public Health (AREA)
- Life Sciences & Earth Sciences (AREA)
- Engineering & Computer Science (AREA)
- Medical Informatics (AREA)
- General Health & Medical Sciences (AREA)
- Biomedical Technology (AREA)
- Veterinary Medicine (AREA)
- Animal Behavior & Ethology (AREA)
- Pathology (AREA)
- Surgery (AREA)
- Data Mining & Analysis (AREA)
- Heart & Thoracic Surgery (AREA)
- Biophysics (AREA)
- Physics & Mathematics (AREA)
- Neurosurgery (AREA)
- Neurology (AREA)
- Molecular Biology (AREA)
- Databases & Information Systems (AREA)
- Epidemiology (AREA)
- Primary Health Care (AREA)
- Nuclear Medicine, Radiotherapy & Molecular Imaging (AREA)
- Radiology & Medical Imaging (AREA)
- Measurement And Recording Of Electrical Phenomena And Electrical Characteristics Of The Living Body (AREA)
Abstract
本发明公开一种自适应闭环深部脑刺激方法、装置及电子设备，其中方法包括：通过粒子群优化算法进行参数搜索获取目标比例-微分-积分控制器参数，采用与目标比例-微分-积分控制器参数对应的刺激参数进行深部脑刺激，进行深部脑刺激的同时对目标比例-微分-积分控制器参数进行在线校正，该方法可实现针对不同患者自动计算控制器增益且随着患者状态的改变自动校准控制器增益，以实现多状态、长时程条件下的个性化神经调控。
Description
本发明涉及医疗电子系统领域，尤其涉及一种自适应闭环深部脑刺激方法、装置及电子设备。
大脑内部神经活动的异常会导致多种神经和精神疾病，深部脑刺激技术(deep brain stimulation,DBS)是目前临床上可行的治疗手段，尤其是针对药物治疗无效的患者。
目前临床上使用持续性开环刺激模式，由医生根据患者的情况调节刺激参数并固定，直到下次复诊再次重新调节。刺激参数包括刺激幅值、频率以及脉宽。在开环刺激模式中，刺激参数无法根据患者瞬时或长期的状态变化进行适当的调整，长期持续性的刺激也有可能带来诸如语言障碍、认知功能障碍等许多副作用。
随着闭环DBS技术的成熟，闭环DBS系统通过施加可根据控制目标进行调节的电刺激信号控制大脑的神经活动。如何通过检测病理性神经活动，形成自适应闭环DBS成为了提高DBS临床治疗效果需要解决的重要问题。
发明内容
本发明的目的在于提供一种自适应闭环深部脑刺激方法、装置及电子设备，其能有效解决多状态、长时程条件下的个性化神经调控问题。
为实现上述发明目的，本发明提出了如下技术方案：
第一方面，一种自适应闭环深部脑刺激方法，所述方法包括：
通过粒子群优化算法进行参数搜索获取目标比例-微分-积分控制器参数以确定所述目标比例-微分-积分控制器；
采用通过所述目标比例-微分-积分控制器获得的刺激参数进行深部脑刺激；
进行深部脑刺激的同时对所述目标比例-微分-积分控制器参数进行在线校正。
在一种较佳的实施方式中，所述通过粒子群优化算法进行参数搜索获取目标比例-微分-积分控制器参数以确定所述目标比例-微分-积分控制器，包括：
初始化粒子群优化算法参数；
迭代计算粒子群中任一粒子在任一第一窗长内的当前适应度并更新全局最优适应度；
判断当前次迭代后的粒子群是否符合终止迭代条件；
若是，则将当前次迭代后的全局最优适应度对应的粒子的位置坐标作为目标比例-微分-积分控制器参数以确定所述目标比例-微分-积分控制器。
在一种较佳的实施方式中，所述判断当前次迭代后的粒子群是否符合终止迭代条件，包括：
判断当前迭代次数是否达到预设迭代次数；或，
判断当前次迭代后所述粒子群的平均适应度与当前更新后的所述全局最优适应度是否相等；或，
判断当前次迭代后的所述全局最优适应度与预设目标适应度是否相同。
在一种较佳的实施方式中，所述迭代计算粒子群中任一粒子在任一第一窗长内的当前适应度并更新全局最优适应度，包括：
迭代计算粒子群中任一粒子在任一第一窗长内的当前适应度；
当任一粒子在当前第一窗长的当前适应度小于其在先的任一适应度，则将相应粒子的所述当前适应度作为相应粒子的个体最优适应度；
当任一粒子的个体最优适应度小于所述粒子群中其余粒子的个体最优适应度时，则将相应粒子的所述个体最优适应度作为全局最优适应度。
在一种较佳的实施方式中，所述迭代计算粒子群中任一粒子在任一第一窗长内的当前适应度，包括：
以任一粒子当前次迭代中的当前位置坐标作为比例-微分-积分控制器参数在当前第一窗长内进行当前次深部脑刺激；
采集所述当前次深部脑刺激时的当前神经活动信号；
基于所述当前神经活动信号与预设目标信号获得所述当前适应度。
在一种较佳的实施方式中，所述采用通过所述目标比例-微分-积分控制器获得的刺激参数进行深部脑刺激，包括：
基于所述目标比例-微分-积分控制器参数获得目标刺激参数；
基于所述目标刺激参数形成刺激脉冲进行深部脑刺激。
在一种较佳的实施方式中，所述进行深部脑刺激的同时对所述目标比例-微分-积分控制器参数进行在线校正，包括：
对深部脑刺激过程中的神经活动信号进行监测并判断是否需要调整所述目标比例-微分-积分控制器参数；
若是，则再次通过粒子群优化算法进行参数搜索以更新所述目标比例-微分-积分控制器参数。
在一种较佳的实施方式中，所述对深部脑刺激过程中的神经活动信号进行监测并判断是否需要调整所述目标比例-微分-积分控制器参数，包括：
基于任一第二窗长中深部脑刺激时的神经活动信号与预设目标信号获得相应的稳态误差；
当所述稳态误差超过预设稳态误差阈值的连续第二窗长数量达到预设窗口数，则判断需要调整所述目标比例-微分-积分控制器参数。
在一种较佳的实施方式中，所述方法还包括：获取自适应闭环深部脑刺激对应的目标脑刺激响应，包括：
获取目标时序刺激输入；
基于所述目标时序刺激输入，通过预先构建的脑刺激响应模型生成目标脑刺激响应。
在一种较佳的实施方式中，所述方法还包括：
预先构建所述脑刺激响应模型，构建方法包括：
获取训练样本集，所述训练样本集包括对大脑的至少一组时序刺激输入及对应于每一组所述时序刺激输入的时序真实响应；
以所述时序刺激输入及所述时序真实响应为输入，基于生成对抗网络训练获得所述脑刺激响应模型。
在一种较佳的实施方式中，所述对抗网络包括生成网络及对抗网络；
所述以所述时序刺激输入及所述时序真实响应为输入，基于生成对抗网络训练获得所述脑刺激响应模型，包括：
以所述时序刺激输入为所述生成网络的输入以获得相应的时序生成响应；
以与所述时序刺激输入相应的所述时序真实响应以及所述时序生成响应为所述对抗网络的输入以获得相应的判断结果；
当所述判断结果符合预设条件时，则停止训练并将与所述生成网络对应的模型作为所述脑刺激响应模型。
在一种较佳的实施方式中，所述当所述判断结果符合预设条件时，则停止训练并将与所述生成网络对应的模型作为所述脑刺激响应模型，包括：
当所述判断结果为所述时序生成响应与所述时序真实响应相同，则停止训练并将与所述生成网络对应的模型作为所述脑刺激响应模型；或，
当所述判断结果为所述时序生成响应与所述时序真实响应不同，则继续训练至所述判断结果符合预设阈值时，停止训练并将与所述生成网络对应的模型作为所述脑刺激响应模型。
在一种较佳的实施方式中，所述以所述时序刺激输入及所述时序真实响应为输入，基于生成对抗网络训练获得所述脑刺激响应模型，还包括：
训练中采用反向传播算法更新所述生成网络及对抗网络的权值及偏置。
在一种较佳的实施方式中，在训练获得所述脑刺激响应模型后，所述方法还包括：基于预先获取的测试样本集对所述脑刺激响应模型进行模型评价，包括：
将所述测试样本集中的任一时序刺激输入输入所述脑刺激响应模型并获得相应的测试结果；
计算所述测试样本集中与所述时序刺激输入对应的时序真实响应与相应测试结果之间的皮尔森相关系数，当所述皮尔森相关系数符合预设阈值时，测试通过。
在一种较佳的实施方式中，所述时序刺激输入包括时序刺激幅度及时序刺激频率，所述时序真实响应包括采集到的与所述时序刺激幅度及时序刺激频率对应的时序真实局部场电位信号，所述时序生成响应包括基于所述刺激幅度及刺激频率通过所述生成网络生成的时序生成局部场电位信号。
在一种较佳的实施方式中，在所述以所述时序刺激输入及所述时序真实响应为输入，基于生成对抗网络训练获得所述脑刺激响应模型之前，所述方法还包括对所采集的时序真实响应进行预处理，包括：
对所采集的所述时序真实局部场电位信号去除刺激伪迹；
对去除刺激伪迹的所述时序真实局部场电位信号进行降采样；
对降采样后的所述时序真实局部场电位信号进行滤波处理；
基于滤波处理后的所述时序真实局部场电位信号计算相应的功率时间序列获得预处理后的时序真实局部场电位信号。
在一种较佳的实施方式中，所述方法还包括：基于脑部电信号构建疼痛状态预测模型，包括：
1)对脑部电信号进行预处理，去除质量不佳的信号和信号中的噪声；
2)从时间域、频率域和小波域三个维度提取特征对脑部电活动进行表征；
3)按照与疼痛状态的相关性，在时间域和小波域上筛选脑部电活动特征；在频率域上使用主成分分析法PCA根据贡献率获得表征各特征组的关键成分进而筛选脑部电活动特征；
4)将从时间域、频率域和小波域三个维度筛选出的特征作为自变量，疼痛缓解程度作为因变量，通过回归分析分别建立状态预测模型；
5)将不同维度上的状态预测结果作为自变量，患者的临床主观评估结果作为因变量，利用多元回归分析建立整合性的疼痛状态预测模型；其中：
步骤2)中，时间域提取的特征包括信号幅值的平均值、标准差以及信息熵；在提取特征前对信号进行标准化处理，具体方法为将各采样点的信号值除以幅值的最大值；小波域提取的特征为delta、theta、alpha、low-beta、high-beta、low-gamma和high-gamma频段的同步化状态存在时间的百分比和这7个频率段两两组合得到的21种组合状态中的各频率段的同步化水平的二值化编码组成的4种状态00、01、10、11出现时间占总时间的百分比。
在一种较佳的实施方式中，步骤1)中，对脑部电信号进行预处理包括：
去50Hz工频干扰和基线漂移；
对去除噪声的信号进行归一化处理。
在一种较佳的实施方式中，步骤2)中，频率域提取的特征为傅里叶变换之后功率谱密度在不同频率段上积分的功率值以及不同频带间功率的比值；在提取特征前对信号进行标准化处理，具体方法为将每个频率点上的功率谱密度值除以2-90Hz频率段功率谱密度的积分。
在一种较佳的实施方式中，步骤3)中，在时间域和小波域上，选择疼 痛状态显著性小于0.05或0.01的特征；在频率域上，根据贡献率最大的1-3个主成分选择特征。
在一种较佳的实施方式中，步骤5)中，疼痛状态预测模型每完成一次预测，就将当前的数据纳入其中用以修正模型参数。
第二方面，提供一种自适应闭环深部脑刺激装置，所述装置包括：
参数搜索模块，用于通过粒子群优化算法进行参数搜索获取目标比例-微分-积分控制器参数以确定所述目标比例-微分-积分控制器；
刺激模块，用于采用通过所述目标比例-微分-积分控制器获得的刺激参数进行深部脑刺激；
校正模块，用于进行深部脑刺激的同时对所述目标比例-微分-积分控制器参数进行在线校正。
在一种优选的实施方式中，所述装置还包括：
第一获取模块，用于获取目标时序刺激输入；
生成模块，用于基于所述目标时序刺激输入，通过预先构建的脑刺激响应模型生成目标脑刺激响应。
第三方面，提供一种电子设备，包括：
一个或多个处理器；以及
与所述一个或多个处理器关联的存储器，所述存储器用于存储程序指令，所述程序指令在被所述一个或多个处理器读取执行时，执行如第一方面任意一项所述的方法。
与现有技术相比，本申请具有如下有益效果：
本申请提供一种自适应闭环深部脑刺激方法、装置及电子设备，其中方法包括：通过粒子群优化算法进行参数搜索获取目标比例-微分-积分控制器参数，采用与目标比例-微分-积分控制器参数对应的刺激参数进行深部脑刺激，进行深部脑刺激的同时对目标比例-微分-积分控制器参数进行在线校正，该方法可实现针对不同患者自动计算控制器增益且随着患者状态的改变自动校准控制器增益，以实现多状态、长时程条件下的个性化神经调控；
进一步的，该申请还包括获取自适应闭环深部脑刺激对应的目标脑刺 激响应，包括获取目标时序刺激输入；基于目标时序刺激输入，通过预先构建的脑刺激响应模型生成目标脑刺激响应；该方法能针对大脑在响应刺激过程中的时变性、非线性及不确定性等特性，真实模拟出大脑在受到深部脑刺激之后的响应，以便提高脑刺激参数的准确性；
进一步的，该申请还包括基于脑部电信号构建疼痛状态预测模型，其能从多维角度全面表征量化脑部电活动，将主观的评估和客观的检测手段结合起来，并且将检测出的多个生物标记进行融合，建立患者状态或变化程度的定量预测模型，该模型能用于实现或者疼痛状态的准确判断；
当然，本申请仅需实现上述任一技术效果即可。
图1是本实施例中自适应闭环深部脑刺激方法的流程图；
图2是本实施例中自适应闭环深部脑刺激方法的又一流程图；
图3为本实施例中参数计算与自适应刺激时序示意图；
图4为本实施例中基于PID控制器的闭环深部脑刺激原理示意图；
图5是本实施例中脑刺激响应模型构建方法的流程图；
图6是本实施例中脑刺激响应模型构建方法的示意图；
图7为本实施例中脑刺激响应方法的流程图；
图8是构建基于脑部电信号疼痛状态(疼痛患者的疼痛缓解程度)预测模型的示例流程示意图。
图9是使用功率谱分析对神经病理性疼痛患者脑部内所记录的电信号(即LFP)的频率点的功率和频率带之间功率的比值的计算结果示例图。
图10是频率域特征与患者疼痛缓解程度进行相关分析得到的相关系数结果图，并且示出了识别出的与疼痛患者程度相关的特征。
图11是在疼痛患者数据中得到的预测模型性能比较及验证结果示例图。
为使本发明的目的、技术方案和优点更加清楚，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显 然，所描述的实施例仅仅是本发明一部分实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都属于本发明保护的范围。
在本发明的描述中，需要理解的是，术语“第一”、“第二”仅用于描述目的，而不能理解为指示或暗示相对重要性或者隐含指明所指示的技术特征的数量。由此，限定有“第一”、“第二”的特征可以明示或者隐含地包括一个或者更多个该特征。在本发明的描述中，除非另有说明，“多个”的含义是两个或两个以上。
大脑具有非线性、时变、非平稳的特性，采用传统的脉冲神经元模型和LSSM模型或很难概括出如何驱动个体的大脑网络动力学，或不符合大脑对刺激的响应是非线性的特性，鉴于上述脑刺激现状，需要寻找一种对应于上述大脑特性的脑刺激响应方法，以真实模拟大脑在受到脑刺激时的响应。以下将结合附图1～11对本实施例的脑刺激响应模型构建方法、响应方法、装置及电子设备作进一步的详细描述。
实施例
如图1～4所示，本实施例提供一种自适应闭环深部脑刺激方法，包括如下步骤：
S100、通过粒子群优化算法(Particle Swarm Optimization,PSO)进行参数搜索获取目标比例-微分-积分控制器参数以确定所述目标比例-微分-积分控制器。
PSO是通过模拟鸟群觅食行为而发展起来的一种基于群体协作的随机搜索算法，可用于解决优化问题。将每一只鸟抽象为参数空间中的“粒子”，食物的位置抽象为符合要求的可行解。所有的粒子都有一个由被优化的函数决定的适应值，每个粒子的速度决定他们进行参数搜索的方向和距离，并在每一代计算后更新。迭代过程中，粒子以全局最优解和个体历史最优解为指导在解空间中搜索。
比例-微分-积分(Proportional-Integral-Differential,PID)控制器包括比例环节、积分环节、微分环节共三个控制环节。比例环节即时成比例地反应控制系统的偏差信号e(t)，偏差一旦产生，控制器立即产生控制作用以减小误差。当偏差e＝0时，控制作用也为0。积分环节基于历史误差产生控制作 用，主要用于消除静差，提高系统的准确度，积分作用的强弱取决于积分时间常数，积分环节增益越大，积分作用越弱，反之则越强。微分环节能反映偏差信号的变化趋势(变化速率)，并能在偏差信号值变得太大之前，在系统中引入一个有效的早期修正信号，从而加快系统的动作速度，减小调节时间。PID控制器具有算法简单、鲁棒性好、可靠性高的特点。特别的，本实施例中具体使用的是增量式PID控制器，属于经典PID的一个变型。
具体地，步骤S100包括：
S110、初始化粒子群优化算法参数。
具体地，对表1所示PSO参数进行初始化。其中，粒子位置x与粒子速度v为随机值；粒子位置维数d对应待确定的控制器增益数量；粒子数N、最大迭代次数G、每次迭代的第一窗长t
1、惯性权重w、加速度常数c
1/c
2、初始值范围b为人为设定值。
表1
| 参数名称 | 符号 |
| 粒子位置坐标 | x |
| 粒子速度 | v |
| 粒子数 | N |
| 粒子位置维数 | d |
| 最大迭代次数 | G |
| 每次迭代时的第一窗长 | t 1 |
| 惯性权重 | w |
| 加速度常数 | c 1,c 2 |
| 初始位置范围 | b＝[b u,b l] |
其中，初始位置范围b的取值范围会影响PSO搜索过程的耗时，如不加限制，由于初始位置是随机的，可能距离满足要求的最优解很远，则搜索过程耗时长。
S120、迭代计算粒子群中任一粒子在任一第一窗长内的当前适应度并更新全局最优适应度。具体的，步骤S12包括：
S121、迭代计算粒子群中任一粒子在任一第一窗长内的当前适应度。
需要说明的是，如图3所示，当基于PID控制器的闭环深部脑刺激时，原理如下：
采用PID控制器以一定的增益进行深部脑刺激并在刺激的同时获取神 经活动信号，信号经处理后相较于预设目标信号获得一个适应度(稳态误差)，用于对PID控制器增益进行调整，从而实现实时闭环控制。
S121具体包括：
S121a、以任一粒子当前次迭代中的当前位置坐标作为PID控制器参数在当前第一窗长内进行当前次深部脑刺激。其中，当前位置坐标即粒子位置向量。
具体的，本实施例中将粒子的当前位置坐标x＝(K
p,K
i,K
d)作为当前第一窗长t
1的PID控制器参数。在确定PID控制器参数后，采用如下公式(1)、(2)计算刺激参数：
Δu(k)＝K
p[e(k)-e(k-1)]+K
ie(k)+K
d[e(k)-2e(k-1)+e(k-2)] (2)
u(k)即为刺激参数。通常，u(k)为幅值、频率、脉宽中的至少一种。本实施例中，频率、脉宽优选为设定值，如130Hz、60μs，故u(k)为幅值。
在通过上述公式(1)(2)获取脑刺激幅值后，即进行深部脑刺激，刺激时长为第一窗长t
1。
S121b、刺激同时，采集当前次深部脑刺激时的当前神经活动信号。
S121c、基于当前神经活动信号与预设目标信号获得当前适应度。
通常的，如图4所示，在实际执行S121a～S121c过程中，由于每次均以每一随机粒子在当前窗长的当前位置坐标作为PID增益进行深部脑刺激，故在每一次迭代过程中实际需要在N个第一窗长t
1中的每个窗长依次进行与每一随机粒子的当前位置坐标作为PID增益进行深部脑刺激。即，在完成粒子1在第一窗长t1中的当前位置坐标作为PID增益进行深部脑刺激并获得粒子1在第一窗长t1中适应度之后，依次持续完成粒子2、3……N在按序排列的第一窗长t
1中的当前位置坐标作为PID增益分别进行深部脑刺激并获得相应的适应度，才完成第一次迭代，并在完成第一次迭代之后持续进行后续迭代。
进一步，在开始PID控制器参数搜索后，以随机的粒子1在第一窗长t
1 中的当前位置坐标作为PID增益进行深部脑刺激并计算获得粒子1在第一窗长t
1中的当前适应度。为了让大脑在一次刺激之后摆脱刺激影响而恢复到未被刺激的状态，以免对下一第一窗长t
1进行以粒子2当前坐标作为PID增益进行深部脑刺激造成不必要的影响，作为一种优选的实施方式，在粒子1对应的第一窗长t
1之后，设置如图4所示的洗脱期t
2。在本实施例的迭代过程中，在每一次第一窗长t
1之后均设置洗脱期t
2，在洗脱期t
2期间对大脑不作任何刺激，以保证粒子群优化过程中的每一深部脑刺激的大脑状态是一致的。
需要说明的是，本实施例中通过适应度函数计算相应的适应度。以控制目标仅为抑制某一频段神经活动的能量为例，当粒子的当前位置坐标x＝(K
p,K
i,K
d)，适应度函数为最小化反馈的该频段能量与预设目标信号差值的平方和，如式(3)所示：
其中，t
s表示每个时间窗长的起始时刻，t
e代表每个时间窗长的结束时刻。其中的e(t)表示在当前窗长中每个反馈时刻的误差。最小化该频段能量与预设目标信号的差值，预设目标信号为任意具有生理意义的神经活动特征，可以为预先设置的经验值，也可以为经过算法获得的个性化数值，本实施例对此不作限制。
当然，上述的控制目标可以由单一目标扩展为多目标，当扩展为多目标调控时，粒子位置向量的维数相应变化。
另外，通过修改适应度函数，也可以达到多重性能指标。例如，若要求闭环深部脑刺激系统在稳态误差尽可能小的情况下，控制输入也尽可能小，即以尽量小的刺激能量完成闭环神经调控，则适应度函数可修改为式(4)：
其中，w
1和w
2分别表示误差和控制量对参数选择的惯性权重。另外， 受限于PID控制器对多变量调控的局限性，本方法中控制目标数量过多时可能受限。
S122、当任一粒子在当前第一窗长的当前适应度小于其在先的任一适应度，则将相应粒子的当前适应度作为其个体最优位置pbest的个体最优适应度，即进行单一粒子本身的个体最优位置pbest及个体最优适应度更新。
S123、当任一粒子的个体最优适应度小于粒子群中其余粒子的个体最优适应度时，则将相应粒子的个体最优适应度作为全局最优适应度，即进行粒子群的全局最优适应度更新，该粒子的粒子位置为全局最优位置gbest。
在更新pbest、gbest之后即可更新粒子的速度和粒子位置，具体如公式(5)、(6)所示：
v
i＝ωv
i+c
1rand(0，1)(pbest
i，d-x
i)+c
2rand(0，1)(gbest
d-x
i) (5)
x
i，d-1＝x
i，d+v
i (6)
其中，i＝1,…,N
i为粒子编号，pbest
i,d为第i个粒子在第d次迭代时的个体最优位置，gbest
d为第d次迭代时的全局最优位置，w为惯性权重，c
1、c
2分别为加速度常数。
S130、在上述更新过程中，判断当前次迭代后的粒子群是否符合终止迭代条件；若是，则将当前次迭代后的全局最优适应度对应的粒子的位置坐标作为目标PID控制器参数，从而确定目标比例-微分-积分控制器。本实施例对于终止迭代条件不作限制，包括但不限于迭代次数、平均适应度大小或全局最优适应度中的任意一项的判断结果。
具体地，步骤S130包括：
S130a、判断当前迭代次数是否达到预设迭代次数；或，
S130b、判断当前次迭代后粒子群的平均适应度与当前更新后的全局最优适应度是否相等；或，
S130c、判断当前次迭代后的全局最优适应度与预设目标适应度是否相同。
上述步骤S130a～S130c择一执行，且当迭代过程中满足任一条件即可触发下一步骤S2。当然，在不满足上述任一条件时则持续进行迭代计算。
经上述步骤S13判断后即可根据全局最优位置gbest确定相应粒子位置 坐标x＝(K
p,K
i,K
d)，即确定步骤S2中进行深部脑刺激的目标PID控制器参数。
S200、采用通过目标PID控制器获得的刺激参数进行深部脑刺激，具体包括：
S210、基于目标PID控制器参数获得目标刺激参数。
具体地，在该步骤中，通过公式(1)、(2)计算目标刺激参数，且公式中的K
p,K
i,K
d为步骤S1中获得的全局最优位置gbest确定的相应粒子位置坐标x＝(K
p,K
i,K
d)。
与步骤S121相似的，本实施例中优选频率、脉宽优选为设定值，故目标刺激参数u(k)指幅值。
S220、基于目标刺激参数形成刺激脉冲进行深部脑刺激。
在步骤S2进行深部脑刺激的同时，继续参照图1、2、4所述，该自适应闭环深部脑刺激方法还包括：
S300、进行深部脑刺激的同时对目标PID控制器参数进行在线校正。步骤S30具体包括：
S310、对深部脑刺激过程中的神经活动信号进行监测并判断是否需要调整目标PID控制器参数；包括：
S311、基于任一第二窗长t
3中深部脑刺激完成时的神经活动信号与预设目标信号获得相应的稳态误差；
S312、当稳态误差超过预设稳态误差阈值的连续第二窗长t
3数量达到预设窗口数，则判断需要调整PID控制器参数。
具体的，稳态误差即为前述的适应度。由于被刺激对象的生理状态处于实时变化中，特别是当被刺激对象病情改变、服药、运动等状态改变时，在同一刺激参数下的适应度将明显增大，若持续刺激将对刺激对象造成不良影响。
S320、若是，则再次重复执行通过粒子群优化算法进行参数搜索以更新目标PID控制器参数，并以更新后的刺激参数进行深部脑刺激，并持续在线调整。在持续在线调整及刺激的具体方法可参照上述步骤S1、S2的描述，此处不再赘述。
若否，则继续执行步骤S200。
因此，本实施例提供的自适应闭环深部脑刺激方法可实现针对不同患 者自动计算控制器增益且随着患者状态的改变自动校准控制器增益，以实现多状态、长时程条件下的个性化神经调控。
在此基础上，本实施例中的自适应闭环深部脑刺激方法还包括S400、获取自适应闭环深部脑刺激对应的目标脑刺激响应，以模拟出真实大脑对脑刺激的响应。具体的，如图5、6所示，步骤S400包括如下步骤：
S410、获取训练样本集，所述训练样本集包括对大脑的至少一组时序刺激输入x
k及对应于每一组时序刺激输入x
k的时序真实响应y
k。
具体地，对大脑交付多种模式的刺激，并记录下刺激输入，同时采集大脑的真实响应输出。通常，刺激输入包括刺激幅度U、刺激频率f，真实响应输出为局部场电位(local field power,LFP)信号y。为提高脑刺激过程中刺激输入与响应输出数据的对应性，本实施例中对训练样本集中的数据均赋予时序性，即刺激输入为时序刺激输入x
k，输出的真实响应也为时序真实响应y
k。
进一步，在脑刺激响应模型构建过程中，时序刺激输入x
k包括时序刺激幅度及时序刺激频率，时序真实响应y
k包括采集到的与时序刺激幅度及时序刺激频率对应的时序真实局部场电位信号。
为提高数据处理精度，通常需要对采集的时序真实响应y
k进行数据预处理，因此在步骤S420之前，还包括：S40、对所采集的时序真实响应y
k进行预处理，具体包括：
S401、对所采集的时序真实局部场电位信号去除刺激伪迹。
刺激过程中记录的宽带原始信号含有刺激伪迹，本实施例中，我们采用模板法来除去刺激伪迹。模板法为脑刺激领域常用的刺激伪迹去除方法，本实施例对此不作进一步描述。
S402、对去除刺激伪迹的时序真实局部场电位信号进行降采样。
示例性的，本实施例采用截止频率为100Hz的抗混叠滤波器将时序真实局部场电位信号降采样至200Hz。
S403、对降采样后的时序真实局部场电位信号进行滤波处理。
示例性的，本实施例采用使用通带截止频率1Hz，阻带截止频率0.5Hz的等纹波有限脉冲响应(FIR)滤波器去除漂移，使用阻带截止频率59Hz和61Hz以及通带截止频率58Hz和62Hz的带阻等纹波FIR滤波器去除60Hz 的线噪声，使用阻带截止频率49Hz和51Hz以及通带截止频率48Hz和52Hz的带阻等纹波FIR滤波器消除在刺激频率上任何可能残留的刺激伪迹。
S404、基于滤波处理后的时序真实局部场电位信号计算相应的功率时间序列获得预处理后的时序真实局部场电位信号。
具体的，根据预设时间窗口长度T
w对每个LFP通道按序划分多个时间窗口并计算每一时间窗口内的LFP的平均功率，由此得到了LFP的功率时间序列y
k，即时序真实局部场电位信号。
时序刺激输入x
k可以为预设经验值。示例性的，时序刺激幅度可选值为0V(无刺激)、1.5V或3V，刺激频率f可选值为0Hz(无刺激)、65Hz或130Hz。因此，训练样本集为
其中x
k＝(U
k,f
k)，N≥2。
在确定训练样本集并进行相应的数据预处理之后，进行步骤S2中的模型训练，具体训练如下所述。
S420、以时序刺激输入x
k及时序真实响应y
k为输入，基于GAN训练获得脑刺激响应模型。当然，在开始模型训练之前，该方法还包括：初始化GAN的权重。
需要说明的是，生成对抗网络(generative adversarial network，GAN)是一种深度学习模型，其包含依次设置的生成网络、判别网络两个网络模型。其中，生成网路的任务是生成看起来自然真实的、和原始数据相似的实例。判别网路的任务是判断给定的实例看起来是自然真实的还是人为伪造的。生成对抗网络通过生成响应的非线性生成器与判别响应真伪的判别器相互竞争的方式训练，能够模拟出真实大脑对脑刺激的响应。
结合图6所示，生成网络的输入为时序刺激输入x
k＝{U
k,f
k}，输出为时序生成响应
时序生成响应
包括基于刺激幅度及刺激频率通过生成网络生成的时序生成局部场电位信号。判别网络的输入为由生成网络产生的时序生成响应
和大脑对刺激的时序真实响应y
k。
因此，步骤S420具体包括：
S421、以时序刺激输入x
k为生成网络的输入以获得相应的时序生成响应
S422、以与时序刺激输入x
k相应的时序真实响应y
k以及时序生成响应
为对抗网络的输入以获得相应的判断结果。
S423、当判断结果符合预设条件时，则停止训练并将与生成网络对应的模型作为脑刺激响应模型。
进一步，步骤S423包括：
S423a、当判断结果为时序生成响应
与时序真实响应y
k相同，则停止训练并将与生成网络对应的模型作为脑刺激响应模型；或，
S423b、当判断结果为时序生成响应
与时序真实响应y
k不同，则继续训练至判断结果符合预设阈值时，停止训练并将与生成网络对应的模型作为脑刺激响应模型。该步骤中用于当时序生成响应
与时序真实响应y
k不相同时，评价两者的接近程度。其中，预设阈值优选0.5，即当判断结果越接近0.5，时序生成响应
与时序真实响应y
k越相近，真实性越高。
上述S423a和S423b择一执行，且当满足S423a时，不再执行S423b。
S424、训练中采用反向传播算法更新生成网络及对抗网络的权值及偏置。
在具体训练中，GAN的目标函数如下式(7)所示：
其中，G(x)为生成网络，D(x)为对抗网络。
判别网络更新梯度如下式(8)所示：
生成网络更新梯度如下式(9)所示：
进一步，在训练获得脑刺激响应模型后，方法还包括：S430、基于预先获取的测试样本集对脑刺激响应模型进行模型评价，具体包括：
S431、将测试样本集中的任一时序刺激输入x
k输入脑刺激响应模型并获得相应的测试结果。
与训练样本集类似的，测试样本集包括时序刺激输入x
k及对应的时序真实响应y
k，时序刺激输入x
k包括时序刺激幅度及时序刺激频率。在执行步骤S432之前，作为一种优选的实施方式，需要对时序真实响应y
k进行数 据预处理，预处理过程参照步骤S0中的描述，此处将不作赘述。
S432、计算测试样本集中与时序刺激输入x
k对应的时序真实响应y
k与相应测试结果之间的皮尔森相关系数。当皮尔森相关系数符合预设阈值时，测试通过。
需要说明的是，皮尔森相关系数(Pearson’s correlation coefficient,CC)用于度量两个变量之间的相关程度，其值介于-1与1之间。本实施例中，CC是基于GAN的时序生成响应
与时序真实响应y
k之间线性相关程度的量，数值越大，说明相关程度越高，脑刺激响应模型的精确度越高。本实施例中的CC表达式如下式(10)所示：
其中，Cov()和Var()分别表示时间序列的协方差和方差。
在训练获得高精度的脑刺激响应模型后，步骤S400还包括S440获取自适应闭环深部脑刺激对应的目标脑刺激响应。如图7所示，该脑刺激响应方法包括：
S441、获取目标时序刺激输入；
S442、基于所述目标时序刺激输入，通过如上述步骤S400的脑刺激响应模型生成目标脑刺激响应。
同样的，上述目标时序刺激输入包括目标刺激幅度U、目标刺激频率f。
通过该脑刺激响应模型所获得的与目标时序刺激输入对应的目标脑刺激响应能有效模拟真实的LFP信号。
因此，本实施例提供的脑刺激响应模型构建方法基于生成对抗网络进行建模获所获得的脑刺激响应模型，基于深度学习强大的学习能力，该脑刺激响应模型能针对大脑在响应刺激过程中的时变性、非线性及不确定性等特性，真实模拟出大脑在受到刺激之后的刺激响应。
在此基础上，本实施例提供的自适应闭环深部脑刺激方法还包括基于脑部电信号构建疼痛状态预测模型，如图8所示，包括：
1)记录患者一定时长的脑部电信号，信号可以是头皮电极记录到的头皮脑电(EEG)，由置于颅骨下接触大脑皮层的皮层脑电图(ECoG)信号或者是由 植入到与疼痛相关的深部核团的场电位信号(LFPs)。记录时间不少于180s。同时还需要记录由临床医护人员对疼痛患者进行状态的评估结果。常见的评估量表为疼痛模拟视觉量表(VAS)。
2)对记录到的脑部电活动进行预处理，首先在原始数据中截取一段较平稳没有奇异值的50s数据。然后利用切比雪夫II型滤波器进行2-90Hz的通带滤波，并使用该滤波器去掉50Hz的工频干扰。最后对信号进行标准化处理，对不同维度特征的提取有所不同，具体的请参见以下内容。
3)从时间域、频率域和小波域三个维度提取特征对脑部活动进行表征。
其中时间域提取的特征包括信号幅值的平均值、标准差以及信息熵。在提取特征前对信号进行标准化处理。具体方法为将各采样点的信号值除以幅值的最大值。
其中频率域特征为傅里叶变换之后对功率谱密度在不同频率段上积分的功率值以及不同频带间功率的比值。功率谱密度在某一频率区间的积分表征该频率段信号的活动程度，不同频率段之间活动程度的比值也作为脑部活动的特征。并且为了消除患者信号间的个体差异性，对每位患者的功率谱进行标准化处理，处理方法是将每个频率点上的功率谱密度值除以2-90Hz频率段功率谱密度的积分；不同频带间功率的比值指计算不同频带信号的幅度值之间的比值，为了寻找恰当的频率组合，本实施例中，采用遍历频率组合的方法，考虑到普遍采用的节律频带一般为4Hz或者其倍数，因此分析采用4Hz的频带0.5Hz的步长来遍历频率组合。
其中小波域特征采用申请人专利201610487800.X中设计的方法共对7个频段的信号进行了同步化状态判别。这些频段分别是delta(3-6Hz)，theta(6-9Hz)，alpha(9-12Hz)，low-beta(12-24Hz),high-beta(24-36Hz),low-gamma(36-60Hz),high-gamma(60-90Hz)。并将7个频率段信号两两组合得到20中组合状态，每种组合共有四种状态，分别为00,01,10,11。分别计算单个频率过度同步化状态占总时间的百分比和组合状态四种状态分别占总时间的百分比。
4)将三个维度得到的所有特征与临床评估的结果进行相关性分析，示例中评估结果与术后12月疼痛减缓程度。基于显著性统计学分析p<0.05提取出在各维度上与患者状态相关的脑部电活动特征。时间域和小波域各特征与 术后疼痛评分和术后减缓程度的相关性分析结果见表1-3。在时间域仅幅值的平均值与术后12月疼痛程度相关，而幅值的信息熵与术后12月疼痛缓解程度相关(表1)。在小波域单个节律只是highbeta的同步化时间与术后12月疼痛减缓程度显著相关。同时由两个节律组成的组合状态中有多个状态与术后12月疼痛减缓程度相关。
如下，表1是示出了时间域特征与神经病理性疼痛患者术后12月疼痛程度和疼痛减缓程度的相关性分析结果。表2是示出了小波域特征7个频段过度同步化状态存在时间百分比与神经病理性疼痛患者术后12月疼痛程度和疼痛减缓程度的相关性分析结果。表3是示出了小波域特征7个频段两两组合的同步化状态存在时间百分比与神经病理性疼痛患者术后12月疼痛减缓程度的显著相关的部分结果。
表2
表3
表4
5)由于频率域提取特征较多，使用PCA分别对频率段功率和频率间功率比值对特征向量进行分析，根据贡献率选择合适的主成分作为表征各特征组的关键成分。图9是使用功率谱分析对神经病理性疼痛患者脑部内所记录的电信号(即LFP)的频率点的功率和频率带间功率的比值的计算结果示例图。
图10是频率域特征与患者疼痛缓解程度进行相关分析得到的相关系数结果图，并且示出了识别出的与疼痛患者程度相关的特征。图10左侧的图中粗实线所在频率的功率谱特征与术后12月疼痛减缓程度显著相关(P<0.01)，而图10右侧图中用方框框起来部分的频带之间功率谱密度的比值特征与术后12月疼痛减缓程度显著相关(P<0.01)。在建模中利用主成分分析的方法仅选择图9和图10中的第一主成分。
6)将三个维度选择的特征作为自变量，术后12月疼痛程度或者缓解程度作为因变量，通过回归分析分别构建预测模型。将预测模型的结果作为自变量，术后12月疼痛程度或者疼痛缓解程度作为因变量，通过多元线性回归分析建立整合三个维度特征的疼痛状态预测模型。图11显示了具体的预测结果。结果表明利用整合时间域、频率域和小波域三个域的特征进行预测效果要优于仅用单个维度的特征。对于术后12月疼痛程度的整合性预测可以达到75％，而对术后12月疼痛减缓程度的整合性预测准确性达到了83％。
因此，本实施例从多维角度全面表征量化脑部电活动，将主观的评估和客观的检测手段结合起来，并且将多个生物标记从单一的脑部电活动中分离出来又进行融合，实现状态的准确判断，建立患者状态或变化程度的定量预测模型。通过本发明的模型可对疼痛状态进行有效预测。
对应于该自适应闭环深部脑刺激方法，本实施例进一步提供一种自适应闭环深部脑刺激装置，包括：
参数搜索模块，用于通过粒子群优化算法进行参数搜索获取目标比例-微分-积分控制器参数以确定目标比例-微分-积分控制器。进一步，参数搜 索模块包括：
初始化单元，用于初始化粒子群优化算法参数。
更新单元，用于迭代计算粒子群中任一粒子在任一第一窗长内的当前适应度并更新全局最优适应度。
具体地，更新单元包括：
第一计算子单元，用于迭代计算粒子群中任一粒子在任一第一窗长内的当前适应度。
具体地，第一计算子单元用于：
以任一粒子当前次迭代中的当前位置坐标作为比例-微分-积分控制器参数在当前第一窗长内进行当前次深部脑刺激；
采集完成所述当前次深部脑刺激时的当前神经活动信号；
基于所述当前神经活动信号与预设目标信号获得所述当前适应度。
第一更新子单元，用于当任一粒子在当前第一窗长的当前适应度小于其在先的任一适应度，则将相应粒子的所述当前适应度作为相应粒子的个体最优适应度。
第二更新子单元，用于当任一粒子的个体最优适应度小于所述粒子群中其余粒子的个体最优适应度时，则将相应粒子的所述个体最优适应度作为全局最优适应度。
第一判断单元，用于判断当前次迭代后的粒子群是否符合终止迭代条件。具体地，第一判断单元用于：
判断当前迭代次数是否达到预设迭代次数；或，
判断当前次迭代后所述粒子群的平均适应度与当前更新后的所述全局最优适应度是否相等；或，
判断当前次迭代后的所述全局最优适应度与预设目标适应度是否相同。
第一更新单元，用于当判断结果为是，则将当前次迭代后的全局最优适应度对应的粒子的位置坐标作为目标比例-微分-积分控制器参数以确定所述目标比例-微分-积分控制器。
刺激模块，用于采用通过目标比例-微分-积分控制器获得的刺激参数进行深部脑刺激。具体地，刺激模块包括：
计算单元，用于基于所述目标比例-微分-积分控制器参数获得目标刺激 参数；
刺激单元，用于基于所述目标刺激参数形成刺激脉冲进行深部脑刺激。
所述装置还包括校正模块，用于进行深部脑刺激的同时对目标比例-微分-积分控制器参数进行在线校正。具体地，校正模块包括：
第二判断单元，用于对深部脑刺激过程中的神经活动信号进行监测并判断是否需要调整所述目标比例-微分-积分控制器参数；
若是，则再次通过粒子群优化算法进行参数搜索以更新所述目标比例-微分-积分控制器参数。
第二判断单元包括：
第二计算子单元，用于基于任一第二窗长中深部脑刺激时的神经活动信号与预设目标信号获得相应的稳态误差；
判断子单元，用于当所述稳态误差超过预设稳态误差阈值的连续第二窗长数量达到预设窗口数，则判断需要调整所述目标比例-微分-积分控制器参数。
进一步，该自适应闭环深部脑刺激装置还包括：
第一获取模块，用于获取训练样本集，所述训练样本集包括对大脑的至少一组时序刺激输入及对应于每一组所述时序刺激输入的时序真实响应；
第一训练模块，用于以所述时序刺激输入及所述时序真实响应为输入，基于生成对抗网络训练获得所述脑刺激响应模型。
其中，所述生成对抗网络包括生成网络及对抗网络。
第一训练模块包括：
生成网络模型单元，用于以所述时序刺激输入为所述生成网络的输入以获得相应的时序生成响应；
对抗网络模型单元，用于以与所述时序刺激输入相应的所述时序真实响应以及所述时序生成响应为所述对抗网络的输入以获得相应的判断结果；当所述判断结果符合预设条件时，则停止训练并将与所述生成网络对应的模型作为所述脑刺激响应模型。
第二更新单元，用于训练中采用反向传播算法更新所述生成网络及对抗网络的权值及偏置。
进一步，对抗网络模型单元具体用于：
当所述判断结果为所述时序生成响应与所述时序真实响应相同，则停止训练并将与所述生成网络对应的模型作为所述脑刺激响应模型；或，
当所述判断结果为所述时序生成响应与所述时序真实响应不同，则继续训练至所述判断结果符合预设阈值时，停止训练并将与所述生成网络对应的模型作为所述脑刺激响应模型。
评价模块，用于基于预先获取的测试样本集对所述脑刺激响应模型进行模型评价，包括：
测试单元，用于将所述测试样本集中的任一时序刺激输入输入所述脑刺激响应模型并获得相应的测试结果；
处理单元，用于计算所述测试样本集中与所述时序刺激输入对应的时序真实响应与相应测试结果之间的皮尔森相关系数，当所述皮尔森相关系数符合预设阈值时，测试通过。
进一步，所述时序刺激输入包括时序刺激幅度及时序刺激频率，所述时序真实响应包括采集到的与所述时序刺激幅度及时序刺激频率对应的时序真实局部场电位信号，所述时序生成响应包括基于所述刺激幅度及刺激频率通过所述生成网络生成的时序生成局部场电位信号。
进一步，该脑刺激响应模型构建装置还包括预处理模块，用于对所采集的时序真实响应进行预处理。预处理模块包括：
第一预处理单元，用于对所采集的所述时序真实局部场电位信号去除刺激伪迹；
第二预处理单元，对去除刺激伪迹的所述时序真实局部场电位信号进行降采样；
第三预处理单元，对降采样后的所述时序真实局部场电位信号进行滤波处理；
第四预处理单元，基于滤波处理后的所述时序真实局部场电位信号计算相应的功率时间序列获得预处理后的时序真实局部场电位信号。
该装置还包括：
第二获取模块，用于获取目标时序刺激输入；
生成模块，用于基于所述目标时序刺激输入，通过脑刺激响应模型生成目标脑刺激响应。
进一步的，该自适应闭环深部脑刺激装置还包括：
第二训练模块，包括：
第五预处理单元，用于对脑部电信号进行预处理，去除质量不佳的信号和信号中的噪声。预处理包括去50Hz工频干扰和基线漂移的步骤，以及对去除噪声的信号进行归一化处理的步骤。
表征单元，用于从时间域、频率域和小波域三个维度提取特征对脑部电活动进行表征。
具体的，时间域提取的特征包括信号幅值的平均值、标准差以及信息熵；在提取特征前对信号进行标准化处理，具体方法为将各采样点的信号值除以幅值的最大值；小波域提取的特征为delta、theta、alpha、low-beta、high-beta、low-gamma和high-gamma频段的同步化状态存在时间的百分比和这7个频率段两两组合得到的21种组合状态中的各频率段的同步化水平的二值化编码组成的4种状态00、01、10、11出现时间占总时间的百分比。
频率域提取的特征为傅里叶变换之后功率谱密度在不同频率段上积分的功率值以及不同频带间功率的比值；在提取特征前对信号进行标准化处理，具体方法为将每个频率点上的功率谱密度值除以2-90Hz频率段功率谱密度的积分。
筛选单元，用于按照与疼痛状态的相关性，在时间域和小波域上筛选脑部电活动特征；在频率域上使用主成分分析法PCA根据贡献率获得表征各特征组的关键成分进而筛选脑部电活动特征。在时间域和小波域上，选择疼痛状态显著性小于0.05或0.01的特征；在频率域上，根据贡献率最大的1-3个主成分选择特征。
模型构建单元，用于将从时间域、频率域和小波域三个维度筛选出的特征作为自变量，疼痛缓解程度作为因变量，通过回归分析分别建立状态预测模型；将不同维度上的状态预测结果作为自变量，患者的临床主观评估结果作为因变量，利用多元回归分析建立整合性的疼痛状态预测模型。疼痛状态预测模型每完成一次预测，就将当前的数据纳入其中用以修正模型参数。
需要说明的是：上述实施例提供的自适应闭环深部脑刺激装置在触发自适应闭环深部脑刺激业务时，仅以上述各功能模块的划分进行举例说明，实际应用中，可以根据需要而将上述功能分配由不同的功能模块完成，即 将系统的内部结构划分成不同的功能模块，以完成以上描述的全部或者部分功能。另外，上述实施例提供的自适应闭环深部脑刺激装置与自适应闭环深部脑刺激方法的实施例属于同一构思，即该系统是基于该方法的，其具体实现过程详见方法实施例，这里不再赘述。
另外，本实施例还提供一种电子设备，包括：
一个或多个处理器；以及
与所述一个或多个处理器关联的存储器，所述存储器用于存储程序指令，所述程序指令在被所述一个或多个处理器读取执行时，执行前述的自适应闭环深部脑刺激方法。
关于执行程序指令所执行的自适应闭环深部脑刺激方法，具体执行细节及相应的有益效果与前述方法中的描述内容是一致的，此处将不再赘述。
上述所有可选技术方案，可以采用任意结合形成本发明的可选实施例，即可将任意多个实施例进行组合，从而获得应对不同应用场景的需求，均在本申请的保护范围内，在此不再一一赘述。
需要说明的是，以上所述仅为本发明的较佳实施例，并不用以限制本发明，凡在本发明的精神和原则之内，所作的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。
Claims (24)
- 一种自适应闭环深部脑刺激方法，其特征在于，所述方法包括：通过粒子群优化算法进行参数搜索获取目标比例-微分-积分控制器参数以确定所述目标比例-微分-积分控制器；采用通过所述目标比例-微分-积分控制器获得的刺激参数进行深部脑刺激；进行深部脑刺激的同时对所述目标比例-微分-积分控制器参数进行在线校正。
- 如权利要求1所述的方法，其特征在于，所述通过粒子群优化算法进行参数搜索获取目标比例-微分-积分控制器参数以确定所述目标比例-微分-积分控制器，包括：初始化粒子群优化算法参数；迭代计算粒子群中任一粒子在任一第一窗长内的当前适应度并更新全局最优适应度；判断当前次迭代后的粒子群是否符合终止迭代条件；若是，则将当前次迭代后的全局最优适应度对应的粒子的位置坐标作为目标比例-微分-积分控制器参数以确定所述目标比例-微分-积分控制器。
- 如权利要求2所述的方法，其特征在于，所述判断当前次迭代后的粒子群是否符合终止迭代条件，包括：判断当前迭代次数是否达到预设迭代次数；或，判断当前次迭代后所述粒子群的平均适应度与当前更新后的所述全局最优适应度是否相等；或，判断当前次迭代后的所述全局最优适应度与预设目标适应度是否相同。
- 如权利要求2所述的方法，其特征在于，所述迭代计算粒子群中任一粒子在任一第一窗长内的当前适应度并更新全局最优适应度，包括：迭代计算粒子群中任一粒子在任一第一窗长内的当前适应度；当任一粒子在当前第一窗长的当前适应度小于其在先的任一适应度，则将相应粒子的所述当前适应度作为相应粒子的个体最优适应度；当任一粒子的个体最优适应度小于所述粒子群中其余粒子的个体最优 适应度时，则将相应粒子的所述个体最优适应度作为全局最优适应度。
- 如权利要求4所述的方法，其特征在于，所述迭代计算粒子群中任一粒子在任一第一窗长内的当前适应度，包括：以任一粒子当前次迭代中的当前位置坐标作为比例-微分-积分控制器参数在当前第一窗长内进行当前次深部脑刺激；采集所述当前次深部脑刺激时的当前神经活动信号；基于所述当前神经活动信号与预设目标信号获得所述当前适应度。
- 如权利要求1所述的方法，其特征在于，所述采用通过所述目标比例-微分-积分控制器获得的刺激参数进行深部脑刺激，包括：基于所述目标比例-微分-积分控制器参数获得目标刺激参数；基于所述目标刺激参数形成刺激脉冲进行深部脑刺激。
- 如权利要求1～6任意一项所述的方法，其特征在于，所述进行深部脑刺激的同时对所述目标比例-微分-积分控制器参数进行在线校正，包括：对深部脑刺激过程中的神经活动信号进行监测并判断是否需要调整所述目标比例-微分-积分控制器参数；若是，则再次通过粒子群优化算法进行参数搜索以更新所述目标比例-微分-积分控制器参数。
- 如权利要求7所述的方法，其特征在于，所述对深部脑刺激过程中的神经活动信号进行监测并判断是否需要调整所述目标比例-微分-积分控制器参数，包括：基于任一第二窗长中深部脑刺激时的神经活动信号与预设目标信号获得相应的稳态误差；当所述稳态误差超过预设稳态误差阈值的连续第二窗长数量达到预设窗口数，则判断需要调整所述目标比例-微分-积分控制器参数。
- 如权利要求1所述的方法，其特征在于，所述方法还包括：获取自适应闭环深部脑刺激对应的目标脑刺激响应，包括：获取目标时序刺激输入；基于所述目标时序刺激输入，通过预先构建的脑刺激响应模型生成目标脑刺激响应。
- 如权利要求9所述的方法，其特征在于，所述方法还包括：预先构建所述脑刺激响应模型，构建方法包括：获取训练样本集，所述训练样本集包括对大脑的至少一组时序刺激输入及对应于每一组所述时序刺激输入的时序真实响应；以所述时序刺激输入及所述时序真实响应为输入，基于生成对抗网络训练获得所述脑刺激响应模型。
- 如权利要求10所述的构建方法，其特征在于，所述对抗网络包括生成网络及对抗网络；所述以所述时序刺激输入及所述时序真实响应为输入，基于生成对抗网络训练获得所述脑刺激响应模型，包括：以所述时序刺激输入为所述生成网络的输入以获得相应的时序生成响应；以与所述时序刺激输入相应的所述时序真实响应以及所述时序生成响应为所述对抗网络的输入以获得相应的判断结果；当所述判断结果符合预设条件时，则停止训练并将与所述生成网络对应的模型作为所述脑刺激响应模型。
- 如权利要求11所述的构建方法，其特征在于，所述当所述判断结果符合预设条件时，则停止训练并将与所述生成网络对应的模型作为所述脑刺激响应模型，包括：当所述判断结果为所述时序生成响应与所述时序真实响应相同，则停止训练并将与所述生成网络对应的模型作为所述脑刺激响应模型；或，当所述判断结果为所述时序生成响应与所述时序真实响应不同，则继续训练至所述判断结果符合预设阈值时，停止训练并将与所述生成网络对应的模型作为所述脑刺激响应模型。
- 如权利要求11所述的构建方法，其特征在于，所述以所述时序刺激输入及所述时序真实响应为输入，基于生成对抗网络训练获得所述脑刺激响应模型，还包括：训练中采用反向传播算法更新所述生成网络及对抗网络的权值及偏置。
- 如权利要求11所述的方法，其特征在于，在训练获得所述脑刺激响 应模型后，所述方法还包括：基于预先获取的测试样本集对所述脑刺激响应模型进行模型评价，包括：将所述测试样本集中的任一时序刺激输入输入所述脑刺激响应模型并获得相应的测试结果；计算所述测试样本集中与所述时序刺激输入对应的时序真实响应与相应测试结果之间的皮尔森相关系数，当所述皮尔森相关系数符合预设阈值时，测试通过。
- 如权利要求11～14任意一项所述的构建方法，其特征在于，所述时序刺激输入包括时序刺激幅度及时序刺激频率，所述时序真实响应包括采集到的与所述时序刺激幅度及时序刺激频率对应的时序真实局部场电位信号，所述时序生成响应包括基于所述刺激幅度及刺激频率通过所述生成网络生成的时序生成局部场电位信号。
- 如权利要求15所述的构建方法，其特征在于，在所述以所述时序刺激输入及所述时序真实响应为输入，基于生成对抗网络训练获得所述脑刺激响应模型之前，所述方法还包括对所采集的时序真实响应进行预处理，包括：对所采集的所述时序真实局部场电位信号去除刺激伪迹；对去除刺激伪迹的所述时序真实局部场电位信号进行降采样；对降采样后的所述时序真实局部场电位信号进行滤波处理；基于滤波处理后的所述时序真实局部场电位信号计算相应的功率时间序列获得预处理后的时序真实局部场电位信号。
- 如权利要求1所述的方法，其特征在于，所述方法还包括：基于脑部电信号构建疼痛状态预测模型，包括：1)对脑部电信号进行预处理，去除质量不佳的信号和信号中的噪声；2)从时间域、频率域和小波域三个维度提取特征对脑部电活动进行表征；3)按照与疼痛状态的相关性，在时间域和小波域上筛选脑部电活动特征；在频率域上使用主成分分析法PCA根据贡献率获得表征各特征组的关键成分进而筛选脑部电活动特征；4)将从时间域、频率域和小波域三个维度筛选出的特征作为自变量，疼痛缓解程度作为因变量，通过回归分析分别建立状态预测模型；5)将不同维度上的状态预测结果作为自变量，患者的临床主观评估结果作为因变量，利用多元回归分析建立整合性的疼痛状态预测模型；其中：步骤2)中，时间域提取的特征包括信号幅值的平均值、标准差以及信息熵；在提取特征前对信号进行标准化处理，具体方法为将各采样点的信号值除以幅值的最大值；小波域提取的特征为delta、theta、alpha、low-beta、high-beta、low-gamma和high-gamma频段的同步化状态存在时间的百分比和这7个频率段两两组合得到的21种组合状态中的各频率段的同步化水平的二值化编码组成的4种状态00、01、10、11出现时间占总时间的百分比。
- 根据权利要求17所述的方法，其特征在于，步骤1)中，对脑部电信号进行预处理包括：去50Hz工频干扰和基线漂移；对去除噪声的信号进行归一化处理。
- 根据权利要求17所述的方法，其特征在于，步骤2)中，频率域提取的特征为傅里叶变换之后功率谱密度在不同频率段上积分的功率值以及不同频带间功率的比值；在提取特征前对信号进行标准化处理，具体方法为将每个频率点上的功率谱密度值除以2-90Hz频率段功率谱密度的积分。
- 根据权利要求17所述的方法，其特征在于，步骤3)中，在时间域和小波域上，选择疼痛状态显著性小于0.05或0.01的特征；在频率域上，根据贡献率最大的1-3个主成分选择特征。
- 如权利要求17所述的方法，其特征在于，步骤5)中，疼痛状态预测模型每完成一次预测，就将当前的数据纳入其中用以修正模型参数。
- 一种自适应闭环深部脑刺激装置，其特征在于，所述装置包括：参数搜索模块，用于通过粒子群优化算法进行参数搜索获取目标比例-微分-积分控制器参数以确定所述目标比例-微分-积分控制器；刺激模块，用于采用通过所述目标比例-微分-积分控制器获得的刺激参数进行深部脑刺激；校正模块，用于进行深部脑刺激的同时对所述目标比例-微分-积分控制 器参数进行在线校正。
- 如权利要求22所述的装置，其特征在于，所述装置还包括：第一获取模块，用于获取目标时序刺激输入；生成模块，用于基于所述目标时序刺激输入，通过预先构建的脑刺激响应模型生成目标脑刺激响应。
- 一种电子设备，其特征在于，包括：一个或多个处理器；以及与所述一个或多个处理器关联的存储器，所述存储器用于存储程序指令，所述程序指令在被所述一个或多个处理器读取执行时，执行如权利要求1～21任意一项所述的方法。
Applications Claiming Priority (4)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN202111101877.6 | 2021-09-18 | ||
| CN202111101876.1 | 2021-09-18 | ||
| CN202111101877.6A CN113941090B (zh) | 2021-09-18 | 2021-09-18 | 一种自适应闭环深部脑刺激方法、装置及电子设备 |
| CN202111101876.1A CN113936806B (zh) | 2021-09-18 | 2021-09-18 | 脑刺激响应模型构建方法、响应方法、装置及电子设备 |
Publications (1)
| Publication Number | Publication Date |
|---|---|
| WO2023040521A1 true WO2023040521A1 (zh) | 2023-03-23 |
Family
ID=85602403
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| PCT/CN2022/111364 WO2023040521A1 (zh) | 2021-09-18 | 2022-08-10 | 一种自适应闭环深部脑刺激方法、装置及电子设备 |
Country Status (1)
| Country | Link |
|---|---|
| WO (1) | WO2023040521A1 (zh) |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN118332521A (zh) * | 2024-06-17 | 2024-07-12 | 山东建筑大学 | 一种基于粒子群优化随机森林的地壳形变时序模拟方法 |
| CN119231431A (zh) * | 2024-12-03 | 2024-12-31 | 武汉佳华创新电气有限公司 | 基于多源数据的配电网继电保护整定方法及系统 |
Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101816822A (zh) * | 2010-05-27 | 2010-09-01 | 天津大学 | 功能性电刺激pid参数双源特征融合微粒子群整定方法 |
| WO2016193440A1 (en) * | 2015-06-03 | 2016-12-08 | Cortec Gmbh | System and method for neurostimulation |
| WO2020152359A1 (en) * | 2019-01-24 | 2020-07-30 | Bios Health Limited | System and method for neural control |
| US10856803B1 (en) * | 2017-08-21 | 2020-12-08 | Aqeel, LLC | Method and apparatus for closed-loop brain stimulation |
| CN113936806A (zh) * | 2021-09-18 | 2022-01-14 | 复旦大学 | 脑刺激响应模型构建方法、响应方法、装置及电子设备 |
| CN113941090A (zh) * | 2021-09-18 | 2022-01-18 | 复旦大学 | 一种自适应闭环深部脑刺激方法、装置及电子设备 |
-
2022
- 2022-08-10 WO PCT/CN2022/111364 patent/WO2023040521A1/zh active Application Filing
Patent Citations (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN101816822A (zh) * | 2010-05-27 | 2010-09-01 | 天津大学 | 功能性电刺激pid参数双源特征融合微粒子群整定方法 |
| WO2016193440A1 (en) * | 2015-06-03 | 2016-12-08 | Cortec Gmbh | System and method for neurostimulation |
| US10856803B1 (en) * | 2017-08-21 | 2020-12-08 | Aqeel, LLC | Method and apparatus for closed-loop brain stimulation |
| WO2020152359A1 (en) * | 2019-01-24 | 2020-07-30 | Bios Health Limited | System and method for neural control |
| CN113936806A (zh) * | 2021-09-18 | 2022-01-14 | 复旦大学 | 脑刺激响应模型构建方法、响应方法、装置及电子设备 |
| CN113941090A (zh) * | 2021-09-18 | 2022-01-18 | 复旦大学 | 一种自适应闭环深部脑刺激方法、装置及电子设备 |
Non-Patent Citations (1)
| Title |
|---|
| CHEN, JUNFENG; FAN, XIN-NAN; SU, LI-YUAN: "Self-tuning of PID Parameters Based on Particle Swarm Optimization", JISUANJI-FANGZHEN = COMPUTER SIMULATION, ZHONGGUO HANGTIAN GONGYE ZONGGONGSI, CN, vol. 23, no. 8, 31 August 2006 (2006-08-31), CN , pages 158 - 160, XP009544649, ISSN: 1006-9348 * |
Cited By (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN118332521A (zh) * | 2024-06-17 | 2024-07-12 | 山东建筑大学 | 一种基于粒子群优化随机森林的地壳形变时序模拟方法 |
| CN119231431A (zh) * | 2024-12-03 | 2024-12-31 | 武汉佳华创新电气有限公司 | 基于多源数据的配电网继电保护整定方法及系统 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| US11224747B2 (en) | Closed loop deep brain stimulation systems and methods | |
| WO2023040521A1 (zh) | 一种自适应闭环深部脑刺激方法、装置及电子设备 | |
| US8197395B2 (en) | Pacemaker for treating physiological system dysfunction | |
| US8521294B2 (en) | System and method for dynamically configurable deep brain stimulation | |
| JP2022540861A (ja) | 神経記録の品質をモニタすること | |
| Roberts et al. | Analysis of the sleep EEG using a multilayer network with spatial organisation | |
| JP5282081B2 (ja) | 深部脳刺激療法のアーチファクトをフィルタリングするための方法 | |
| DE102015218948A1 (de) | Vorrichtung und Verfahren zur Bestimmung eines medizinischen Gesundheitsparameters eines Probanden mittels Stimmanalyse | |
| WO2022029445A1 (en) | Deep brain stimulation | |
| CN113112017B (zh) | 基于神经流形的脑电分级与预后fpga解码系统 | |
| Dutta et al. | Oscillation quenching in Stuart-Landau oscillators via dissimilar repulsive coupling | |
| CN113936806A (zh) | 脑刺激响应模型构建方法、响应方法、装置及电子设备 | |
| CN113143275A (zh) | 一种样本与特征质量联合量化评估的脑电疲劳检测方法 | |
| WO2024084391A1 (en) | Method for classifying electroencephalography data with the use of artificial neural networks | |
| Martens et al. | Separating burst from background spikes in multichannel neuronal recordings using return map analysis | |
| CN116528749A (zh) | 用于训练神经疾病检测算法的计算机程序、编程植入式神经刺激装置的方法及其计算机程序 | |
| US12097029B1 (en) | Decoding neuropsychiatric states from multi-site brain network activity | |
| Petkov et al. | Computational model exploration of stimulation based paradigm for detection of epileptic systems | |
| Martin et al. | Differential contributions of subthalamic beta rhythms and neural noise to Parkinson motor symptoms | |
| Sharmila et al. | Integrative Approach for Epileptic Seizure Detection: A Comparative Analysis | |
| CN118522439B (zh) | 一种数字中医理疗方法及系统 | |
| US20220241594A1 (en) | Systems and Methods for Deep Brain Stimulation Using Beta Burst Feedback | |
| US20240342489A1 (en) | Method for iegm-based monitoring of an electrode status of an implantable device | |
| US20230181909A1 (en) | Systems for using local field potential oscillations | |
| US20220241591A1 (en) | Systems and Methods for Deep Brain Stimulation Using Kinematic Feedback for Treatment of Movement Disorders |
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: 22868898 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: 22868898 Country of ref document: EP Kind code of ref document: A1 |
|
| 122 | Ep: pct application non-entry in european phase |
Ref document number: 22868898 Country of ref document: EP Kind code of ref document: A1 |
|
| 32PN | Ep: public notification in the ep bulletin as address of the adressee cannot be established |
Free format text: NOTING OF LOSS OF RIGHTS PURSUANT TO RULE 112(1) EPC (EPO FORM 1205A DATED 24/10/2024) |