WO2021075548A1 - 脳状態推定装置、コンピュータプログラム、脳状態推定方法、脳機能の検査システムおよび方法 - Google Patents

脳状態推定装置、コンピュータプログラム、脳状態推定方法、脳機能の検査システムおよび方法 Download PDF

Info

Publication number
WO2021075548A1
WO2021075548A1 PCT/JP2020/039091 JP2020039091W WO2021075548A1 WO 2021075548 A1 WO2021075548 A1 WO 2021075548A1 JP 2020039091 W JP2020039091 W JP 2020039091W WO 2021075548 A1 WO2021075548 A1 WO 2021075548A1
Authority
WO
WIPO (PCT)
Prior art keywords
brain
data
entropy
subject
activity data
Prior art date
Application number
PCT/JP2020/039091
Other languages
English (en)
French (fr)
Inventor
進 大槻
晃裕 奥野
Original Assignee
株式会社Ｓｐｌｉｎｋ
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 株式会社Ｓｐｌｉｎｋ filed Critical 株式会社Ｓｐｌｉｎｋ
Priority to JP2021552465A priority Critical patent/JPWO2021075548A1/ja
Publication of WO2021075548A1 publication Critical patent/WO2021075548A1/ja

Links

Images

Classifications

    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B10/00—Other methods or instruments for diagnosis, e.g. instruments for taking a cell sample, for biopsy, for vaccination diagnosis; Sex determination; Ovulation-period determination; Throat striking implements
    • A—HUMAN NECESSITIES
    • A61—MEDICAL OR VETERINARY SCIENCE; HYGIENE
    • A61B—DIAGNOSIS; SURGERY; IDENTIFICATION
    • A61B5/00—Measuring for diagnostic purposes; Identification of persons
    • A61B5/05—Detecting, measuring or recording for diagnosis by means of electric currents or magnetic fields; Measuring using microwaves or radio waves 
    • 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
    • G16H10/00—ICT specially adapted for the handling or processing of patient-related medical or healthcare data
    • G16H10/40—ICT specially adapted for the handling or processing of patient-related medical or healthcare data for data related to laboratory analysis, e.g. patient specimen analysis

Definitions

  • the present invention relates to a brain state estimation device, a computer program, a brain state estimation method, a brain function test system and a method.
  • a method of presenting a stimulus to a subject and measuring the reaction of the brain to the stimulus is generally used for a brain function test. It is required to perform appropriate stimulation, measurement, and analysis according to the inspection item by selecting appropriate parameters, and to obtain an appropriate solution for the obtained data.
  • One aspect of the present invention is a system including a function of classifying brain activity data of a subject fluctuating in a space including a time axis by using a time series model in which learning data including a large number of solutions related to a test of brain function is learned.
  • a time series model is a time series neural network including a recurrent neural network (RNN).
  • RNN recurrent neural network
  • One of the other aspects is shown by activity data using a model in which solutions related to a large number of inputs and outputs including brain activity regions and patterns are learned from the subject's brain activity data that fluctuates in space including the time axis. It is a system that includes a function to estimate the spatial distribution including the time axis of the signal source.
  • One of the other aspects is to use a model in which the spatial distribution including the time axis of the signal source estimated from the brain activity data of the subject fluctuating in the space including the time axis and the case are learned to determine the brain state of the subject. It is a system that includes a function to estimate. These models may include time series neural networks.
  • the brain activity data may be resting data and may include brain activity data for stimuli including at least one of visual, auditory, tactile, taste, and olfactory sensations.
  • the present invention includes a program having instructions for operating a computer as the above system.
  • One of the further different aspects of the present invention is to use a time-series model in which the brain activity data of a subject fluctuating in space including the time axis is learned from the learning data including a large number of solutions related to the examination of brain function. It is a method that involves classifying. This method uses a model in which solutions related to a large number of inputs and outputs including brain activity regions and patterns are learned from the subject's brain activity data that fluctuates in space including the time axis, and the signal source indicated by the activity data.
  • It may include estimating the spatial distribution including the time axis of the signal source, and a model that learns the spatial distribution including the time axis of the signal source estimated from the brain activity data of the subject fluctuating in the space including the time axis and the case. It may include estimating the brain state of a subject.
  • Magnetoencephalography is a magnetic field activity in the brain recorded using a magnetoencephalography device (abbreviated as a magnetoencephalogram).
  • EEG electroencephalography
  • MEG magnetic field activity in the brain recorded using a magnetoencephalography device
  • EEG electroencephalography
  • MEG MEG
  • EEG electroencephalogram
  • MEG functional MRI
  • fMRI functional image
  • MEG is currently the most suitable method for obtaining brain activity data, but it may be data obtained by EEG.
  • a neural network is used for MEG analysis, and the state of various brain diseases and the state of brain activity of healthy people are discriminated with high accuracy based on brain activation information, and before analysis such as signal source estimation.
  • a method using a neural network is also created for processing.
  • a process 100 for classifying Source Level data by time series NN may be provided.
  • the time series NN includes a recurrent neural network and the like.
  • a plurality of solution models may be prepared, and a signal source estimation method including a new process 200 for learning a network for reproducing input / output may be provided.
  • a process 300 for discriminating the signal source estimation data using NN may be provided.
  • the object to be differentiated may include brain function and may be differentiated from the presence or absence of symptoms related to brain abnormalities. For example, in dementia, NC / MCI / AD classification may be performed.
  • the brain activity data may be resting data, or may be brain activity data for a stimulus including at least one of visual, auditory, tactile, taste, and olfactory sensations.
  • FIG. 2 shows an experimental example in which resting data of a healthy person was classified into spontaneous eye closure and spontaneous eye opening using a ResNet model.
  • the model that learned the brain activity data related to eye opening and closing it is possible to extract and analyze the data related to eye opening and closing from the brain activity data.
  • FIG. 3 shows, for reference, the radio field intensity of each EEG frequency band regarding the opening and closing of the eyes, and also shows that there is a difference in the radio field strength of each frequency band between healthy subjects and dementia patients. Therefore, by acquiring time-dependently fluctuating brain activity data and using artificial intelligence that has learned knowledge about the appropriate frequency band and knowledge about the source, the soundness of brain function or the function of the brain It is possible to provide a method and a system for accurately determining a change over time for each patient.
  • the space including time is a concept including a space obtained by converting a time domain function into a frequency space by Fourier transform or the like, and the brain obtained by a dynamic image or signal rather than a static image or signal. It relates to systems and methods for handling activity data.
  • Functional connectivity is expressed by changes in the added electrical phenomena (EPSPs: excitatory postsynaptic potentials, IPSPs: inhibitory postsynaptic potentials) caused by networks in the brain made up of neuronal populations.
  • ESPs excitatory postsynaptic potentials
  • IPSPs inhibitory postsynaptic potentials
  • the result of calculating the input / output magnetic field vector strength or potential impedance vector for each region of the anatomical neural projection path is the functional connectivity.
  • Functional connectivity is a type of classification of brain connectivity that expresses statistical relationships based on temporal information on electrical activity projected onto neural populations. Functional connectivity represents the relationships between different regions in the brain as a network by capturing the entire brain from the perspective of brain function.
  • DMC DMC with MNE means a synchronously dispersed neural activity distribution in different brain regions.
  • FIG. 4 is a block diagram showing an example of the configuration of the brain state estimation device 50 of the present embodiment.
  • the brain state estimation device 50 includes a control unit 51 that controls the entire device, a data acquisition unit 52, a storage unit 53, a preprocessing unit 54, an SSP processing unit 55, a signal source position estimation unit 56, a signal distribution calculation unit 57, and a spectrum analysis.
  • a unit 58, an entropy calculation unit 59, a brain state estimation unit 60, a display panel 61, an operation unit 62, and a display processing unit 63 are provided.
  • the control unit 51 can be composed of a CPU (Central Processing Unit), a ROM (Read Only Memory), a RAM (Random Access Memory), and the like.
  • the storage unit 53 can be composed of a semiconductor memory, a hard disk, or the like.
  • MEG data magnetic field data measured by the MEG sensor (detection coil) of the biomagnetic measuring device
  • the brain data means brain activity data, electrooculogram data (EOG), electrocardiogram data (ECG), myoelectric potential data (EMG), reference noise, etc.
  • the brain data includes brain activity data.
  • the brain activity data may include at least one of EOG, ECG, EMG, and reference noise.
  • EOG ECG
  • EMG EMG
  • reference noise can be used in addition to the brain activity data.
  • brain activity data is indispensable when analyzing using the difference in brain state.
  • the brain activity data is not limited to the magnetic field data, but may be data obtained by EEG.
  • the data acquisition unit 52 can acquire the brain data of the subject that fluctuates in the space including the time axis. More specifically, the data acquisition unit 52 can acquire MEG data as brain activity data.
  • the biomagnetic measuring device includes a large number of MEG sensors arranged on the scalp (for example, 160 locations), each MEG sensor corresponds to each channel, and the data acquisition unit 52 supports each of the plurality of channels.
  • MEG data can be acquired.
  • the MEG data can be, for example, data when the subject is in a resting position and the eyes are closed.
  • the data acquisition unit 52 can acquire all or part of data such as electrooculogram (EOG), electrocardiogram (ECG), myoelectric potential (EMG), and reference noise data as data other than MEG data.
  • EEG electroencephalogram
  • EEG electroencephalogram
  • the preprocessing unit 54 performs data removal of defective channels, baseline correction, handpass filter processing, notch filter processing, addition averaging processing, independent component analysis (ICA: Independent Component Analysis), and the like.
  • the measured magnetic field data includes not only magnetic field-derived data but also signals derived from environmental noise, head motion, and biological noise.
  • the preprocessing unit 54 has a function of removing signals other than magnetic field data.
  • Data removal of defective channels removes MEG data of channels that cannot be measured from the data to be analyzed.
  • the determination of the channel that failed to acquire the MEG data is obtained by eigenvalue decomposition for the MEG data group, and the data of the channel showing a predetermined threshold value is removed.
  • the removed channel data can be interpolated (interpolated) with the data of the channel that is not the target of removal.
  • Baseline correction (also referred to as “bias correction”) uses the average value of the time-series data recorded in each channel of the MEG data, and subtracts the average value from the time-series data recorded in each channel. This is because the MEG data is converted into relative data because the absolute value of the signal data fluctuates depending on the temperature, atmospheric pressure, and individual differences.
  • the handpass filter processing uses a bandpass filter having a required frequency band as the bandwidth, and removes signals in an unnecessary frequency band from the MEG data.
  • the required frequency band can be, for example, 0 to 70 Hz.
  • Neurons such as ⁇ wave (8 to 13 Hz), ⁇ wave (14 to 30 Hz), ⁇ wave (4 to 7 Hz), ⁇ wave (0.5 to 3 Hz), and ⁇ wave (26 to 70 Hz) are included in this frequency band. It contains electrophysiological electrical signals of the brain composed of cell populations.
  • the notch filter process uses a notch filter to remove periodic noise (50 Hz or 60 Hz) derived from the power source mixed from the power source of another device or the like used when operating the biomagnetic measuring device.
  • the recorded time series data is divided into, for example, a 5-second interval (trial), and the noise is added and averaged for each divided data. Attenuate the components (epoch data).
  • FIG. 5 is a schematic diagram showing the structure of epoch data.
  • the MEG data of each of the plurality of channels is divided into a plurality of sections (trial1, trial2, ..., Trialn).
  • the time of each section can be, for example, 5 seconds, but is not limited to this. Further, the number n of sections can be appropriately determined.
  • FIG. 6 is a schematic diagram showing an example of addition averaging processing.
  • the signals of each trial from 1 to n are added for each trial.
  • the signals commonly observed in each trial are emphasized by being added.
  • the waveforms randomly observed in each trial are canceled by being added, and the waveform of noise is reduced.
  • Independent component analysis is a method of separating a single signal for each independent signal source from which it is generated.
  • the separated signals include MEG data, environmental noise and EMG data.
  • Independent component analysis uses two pieces of information, a reference channel that measures environmental noise derived from the indoor environment and an EMG channel that measures eye movements and muscle activity, and each is located on the head.
  • the noise component (for example, EMG, noise synchronized with the reference) mixed in the data measured by the MEG sensor is removed. Specifically, a random signal, an artifact associated with head movement, and a signal derived from biological noise are removed from the magnetic field data measured for a single subject, and a plurality of signal components are obtained. Based on the correlation coefficient or threshold, the ICA element of the EMG-EMG, EMG-MEG combination defined as an artifact is removed from each channel.
  • the SSP processing unit 55 reduces the noise effect of the signal source from the MEG data by using the SSP (Signal Space Projection) method.
  • SSP Signal Space Projection
  • the original signal ys (t) at time t is composed of various signal sources (signal sources whose rotation and position are indefinite) and a read field l generated for the signal sources.
  • a certain signal source Sq shows a unique distribution with respect to the channel of MEG, and is represented by the equation (2).
  • q indicates the channel index
  • Q indicates the number of channels
  • lq indicates the signal state of the channel q.
  • the signal source Sq is unknown.
  • the noise P that affects all the influences of the original signal ys (t) and the noise ⁇ can be expressed by the equation (3). When the influence P of noise is minimized, equations (1) to (4) can be defined.
  • the basis vector u consisting of the orthonormal basis of Uq is the maximum likelihood estimation solution of the eigenvector corresponding to the maximum eigenvalue consisting of Q elements of the data covariance matrix.
  • the original signal ys (t) can be obtained by using the singular value decomposition from the covariance matrix of the data in order to obtain the equation (6). This original signal ys (t) is used as sensor space data (Sensor Space Data).
  • the magnetic field data was mixed with independent random number data and EMG + ECG data measured in an independent environment in a magnetic field distribution composed of theoretical voltage values, and the above-mentioned ICA processing and SSP processing were performed.
  • the calculation result of GOF (Goodness of Fit) as an index for evaluating the degree of agreement between the measured value and the theoretical value was as high as 97%.
  • ⁇ means that each sensor is summed.
  • GOF means that when the value is 100%, the two are completely in agreement, and the smaller the value is, the smaller the degree of agreement is. In the case of only ICA treatment, the GOF was as low as 60%.
  • the MEG data can be reconstructed after removing the uncorrelated signals of the time series data of each channel and the uncorrelated signals between the channels.
  • Applicants are related and synchronous in functional connectivity in the brain in healthy subjects, whereas functional connectivity is synchronous in patients with brain diseases such as dementia. I paid attention to the fact that is gradually lost.
  • the correlation component of the signal between the channels is emphasized, and by using the reconstructed MEG data, the healthy subject and the brain disease patient Can be separated (screened).
  • brain diseases include dementia (including AD, DLB, frontal temporal lobe degeneration (FTLD), normal pressure hydrocephalus (NPH), etc.), brain tumors, and psychiatric disorders (also referred to as psychiatric disorders).
  • dementia including AD, DLB, frontal temporal lobe degeneration (FTLD), normal pressure hydrocephalus (NPH), etc.
  • psychiatric disorders also referred to as psychiatric disorders.
  • AD attention deficit / hyperactivity disorder
  • sleep disorder childhood disease
  • ischemic brain disorder includes (including depression, etc.), etc.
  • the signal source position estimation unit 56 estimates the position of the signal source in the brain using the sensor space data obtained by the SSP processing unit 55.
  • the signal source is a current dipole that represents the synchronous activity of a large number of neurons.
  • the magnetic field depends on a function of distance.
  • the measurement data of the MEG sensor is strongly influenced by the distance between the MEG sensor and the head.
  • FIG. 7 is a schematic diagram showing the locality of the brain surface activity distribution accompanying the rotation of the head with respect to the center position of the MEG sensor map.
  • the contour portion indicates the head range that is rotated and expanded according to the head center position and range.
  • the mesh in the contour shows the MNI template rotated and scaled according to the head center position and range.
  • the left and right parts in the contour indicate the area of the shortest distance between the head and the MEG sensor position.
  • the points around the contour represent the position of each MEG sensor.
  • the sensor position In order to reduce the influence of the distance between the MEG sensor and the head, it is necessary to define the sensor position in consideration of the head position. From the spatial positions of the center of the nose, right ear, and left ear recorded in the MEG data, the center position of the head is set, and the affine transformation is performed so that the position of each MEG sensor is set according to the set center position.
  • the head model (MNI template) is rotated and scaled using. Based on the shape of this head model and the sensor distance, the sensor space data obtained by the SSP processing unit 55 is used, and the biot-Savart equation, which is a magnetic field model equation that passes through the medium, is used to signal in the brain. Estimate the source position (dipole estimation by forward problem).
  • the signal source position estimation unit 56 uses the Biot-Savart equation again from the current component at the estimated intracerebral signal source position to form magnetic field data that is specifically attenuated and enhanced by the distance to the sensor located outside the head. Then, the reconstructed sensor space data is used.
  • the signal distribution calculation unit 57 calculates the magnetic field distribution (Source Space Data) using the sensor space data reconstructed by the signal source position estimation unit 56.
  • the magnetic field distribution (SourceSpaceData) is calculated using an adaptive beamformer (SAM).
  • FIG. 8 is a schematic diagram showing a method of calculating the magnetic field distribution composed between the MEG sensor and the current dipole.
  • the magnetic field data of MEG is the magnetic field until it reaches the MEG sensor measured externally when the position and rotation of the dipole Jo with a current component in the brain are given (r, r-r').
  • the distribution A (r) can be calculated using the Biot-Savart equation shown in the figure.
  • SAM estimates the existence of the current dipole Jo using the biosavar equation from the magnetic field data recorded by each MEG sensor, and then uses the mathematical programming method to use an external sensor (Sensor Space Data) or the boundary of the head surface.
  • the magnetic field distribution (Source Space Data) measured on the surface (circular region in the figure) is calculated. By using this method, it is possible to reduce signals other than the magnetic field distribution included in the measured magnetic field data and improve the spatial SNR.
  • FIG. 9 is a schematic diagram showing an example of the result of the magnetic field distribution (Source Space Data) calculated according to the present embodiment.
  • the magnetic field distribution (Source Space Data) obtained by using an adaptive beam former (SAM) improves the spatial SNR and enables accurate evaluation of adjacent neural activity sources.
  • SAM adaptive beam former
  • the spectrum analysis unit 58 can generate information related to the brain by converting the brain data acquired via the data acquisition unit 52. For example, the spectrum analysis unit 58 can use frequency conversion to generate brain-related information, including rhythmic frequency band information.
  • the frequency band information of the rhythm is information that focuses on the fact that the brain function information has a certain periodic structure.
  • the spectrum analysis unit 58 can generate information related to the brain including fractal information by using the transformation by the nonlinear time series analysis.
  • the fractal information is information that focuses on the disorderness contained in the brain function information, and means the degree of breakage of a similar structure of the brain function information.
  • the spectrum analysis unit 58 acquires the magnetic field data (time series data of each MEG channel) obtained by the signal source position estimation unit 56 or the signal distribution calculation unit 57, and performs spectrum analysis of a single channel. , And spectrum analysis between the two channels can be performed.
  • the spectrum analysis can be switched by specifying the setting file including the calculation conditions. In the following, first, a single channel spectrum analysis will be described.
  • the Blackman-Tukey method which obtains the spectrum from the autocorrelation function, is designed to obtain the fluctuation in a finite time domain based on the error theory of the spectrum.
  • a set element consisting of a random sequence in a finite time domain has an individual probability of occurrence in each element.
  • FIG. 10 is a schematic diagram showing an example of the average result in the time domain of two samples having different periodic deviations.
  • the statistic cannot be estimated by using the stochastic process to which a plurality of samples belong, it is necessary to change the time domain to be treated as a sample.
  • a statistic can be shown if all statistics can be estimated from any one sample process x (t) with consistency and impartiality, that is, if they satisfy ergodicity. Therefore, it must be shown that ergodicity is guaranteed by changing the time domain treated as a sample.
  • the time domain T that guarantees the ergodicity of the mean ⁇ and the autocovariance function ⁇ ( ⁇ ) (Autocovariance) is automatically determined.
  • the mean in the time domain is expressed by Eq. (9)
  • the average ergodicity of time series data is expressed by Eq. (10)
  • the zero convergence condition of the variance showing the stationary process is expressed by Eq. (11).
  • C x (t1-t2) represents an autocovariance function in the frequency domains t1 and t2.
  • the power spectral density (PSD) can be estimated.
  • the PSD calculation using the automatically determined time domain T is performed as follows. Assuming that rectangular sampling of time-series data having a certain time domain T is given by the equation (12), the PSD at this time can be obtained by the equation (13).
  • E indicates the time average of x ⁇ (t) x ⁇ (t + ⁇ ) of the time series data x ⁇ in the ⁇ index region normalized by the equation (12).
  • ⁇ indicates the sampling width.
  • FIG. 11 is a schematic diagram showing an example of PSD of the autocovariance function calculated using the time domain T.
  • the horizontal axis represents frequency and the vertical axis represents PSD (power spectral density).
  • PSD power spectral density
  • the solid line shows the theoretical value, and the broken line shows the value using the Blackman-Tukey method in the time domain T.
  • the optimum time domain T can be obtained regardless of the type of window function such as a rectangular window, a Bartlett window, and a Blackman window.
  • FIG. 12 is a schematic diagram showing an example of the difference in PSD between a dementia patient and a healthy person.
  • the horizontal axis represents frequency and the vertical axis represents PSD (power spectral density).
  • the solid line indicates a healthy person (NC), and the broken line indicates an Alzheimer's disease patient (AD).
  • N healthy person
  • AD Alzheimer's disease patient
  • characteristic differences are seen in the ⁇ band and ⁇ band, and the difference in the clinical stage of dementia can be defined by the spectrum. Conversely, it is possible to classify and screen the clinical stage of dementia using spectral density.
  • AD dementia
  • DLB frontal temporal lobe degeneration
  • NPH normal pressure hydrocephalus
  • MCI mild cognitive impairment
  • MCI due to AD prodromal AD
  • presymptomatic stage of Alzheimer's disease / including preclinical AD, etc. multiple sclerosis brain tumors
  • psychiatric disorders also called psychiatric disorders, schizophrenia
  • Parkinson's disease Asperger's syndrome
  • ADHD attention deficit / hyperactivity disorder
  • sleep disorder childhood disease
  • ischemic brain disorder ischemic brain disorder
  • mood disorder Classification discrimination, diagnosis, screening, etc.
  • Equation (19) The inverse Fourier transform of Y'(k), which is a Fourier transform, is the equation (15), and by substituting the equation (15) into the equation (14), the equation (16) is obtained.
  • H (kl) in the equation (16) can be expressed by the equation (17), which is a Fourier transform equation of the taper function h (t).
  • the expected value ⁇ S' ⁇ of the Fourier transform using the taper is expressed by convolution with the Fourier transform of the time series data y (t) and the Fourier transform of the taper function h (t).
  • the expected value ⁇ S' ⁇ of the Fourier transform can be obtained by Eq. (18).
  • f' indicates the sampling frequency.
  • the spectral data using the taper function h (t) is given by Equation (19).
  • the cross spectrum estimation by the combination j between the two channels can be expressed by the equation (20) which extends the equation (16).
  • the cross-spectral estimation is expressed by the equation (20) and the equation (21).
  • the cross spectrum shows a high value when the phases of the data between the two channels for which the cross spectrum is calculated are synchronized.
  • the expected value of the cross spectrum is the same regardless of the Fourier window function (humming window, rectangular window, etc.), and the fluctuation of the cross spectrum estimation error can be made extremely small. ..
  • PSI Phase Slope Index
  • the cross spectrum at this time can be expressed by the equation (24). Assuming that the time at which the phase information is transmitted from the jth position to the kth position is ⁇ , the synchronous phase state of the cross spectrum is proportional to the frequency f of the signal and the time difference ⁇ . Based on this information, Phase Slope ( ⁇ ) is defined by Eq. (25). This can be calculated from the cross spectrum, but can also be described as coherence. Coherence can be described as equation (26) from the cross spectrum ⁇ . Coherence represents the real part of the cross spectrum, which is a complex number, and can be seen as the mutual correlation coefficient of the Fourier frequency components of the two signals. The value shown in equation (26) is between 0 and 1.
  • PSI is determined by equation (27) using coherence. Equation (25) and equation (27) are the same. PSI means synchronization information having a lag in which phase information is transmitted from the jth position to the kth position.
  • the result of spectral analysis using MEG data that is, the frequency band of rhythm (also referred to as “synchronous component, coherence”) coordinated between each region in the brain related to brain activity from the spectral data, in the frequency band of a certain rhythm.
  • the amount of synchronization in the brain (the amount of synchronization phase difference indicated by coherence seen in a certain frequency band) and the synchronization distribution (synchronization in coherence and cross spectrum).
  • the phase difference amount it is possible to calculate the feature amount related to the brain disease (difference in the clinical stage of the brain disease, for example, the entropy that determines the difference between the dementia patient and the healthy person).
  • the entropy calculation unit 59 calculates the entropy based on the information generated by the spectrum analysis unit 58. Specifically, in relation to the two channels, the entropy calculation unit 59 can calculate the entropy that characterizes the synchronous phase difference amount in the brain based on the frequency band information of the rhythm related to the brain activity.
  • Frequency band information includes a rhythmic frequency band and power (also referred to as a "frequency component") in a rhythmic frequency band.
  • the amount of synchronous phase difference refers to the amount of synchronization indicated by coherence seen in a certain frequency band, and the synchronization component in coherence or cross spectrum.
  • Coherence represents the amount of components that are synchronized between the grids, and refers to signal components with the same phase and frequency. Further, in relation to each channel, the entropy calculation unit 59 can calculate the entropy that characterizes the synchronization amount in the brain based on the frequency band information of the rhythm related to the brain activity. In addition, paying attention to the structure of the time-series signal related to the brain activity in relation to each channel, the entropy calculation unit 59 uses the nonlinear time-series analysis to check the disorder of the time-series signal and the degree of disruption of the phase structure. The entropy to characterize can be calculated.
  • the disorder index of the signal and the difference of the time series data of the specified combination are used.
  • Appprocimate entropy (ApEn), which means a measure seen in a certain numerical range
  • Sample Entropy which means a signal randomness index with autocorrelation and does not depend on data size in the extension system of ApEn
  • LZC Koromogorov Lempel-Ziv complexity
  • It means the measurement of time series data using the self-return coefficient.
  • MSE Multisacale entropy
  • RMSE refined MSE
  • Lopez which means the measurement of time series data using self-regression coefficient
  • LMC Ruiz Mancini Calbet complexity
  • the metrics TSE, MF, SE, SSE representing the entropy of a single channel, and the metrics QSE representing the entropy of two channels will be described.
  • the TSE can be calculated by Eq. (28) and is defined by the measure in which a certain spectral density (PSD) is observed in the frequency band of 0.5 Hz to 70 Hz of each channel and the order of the q exponential function. The difference from the most probable function is shown.
  • PSD spectral density
  • the QSE can be calculated by Eq. (29), a measure in which one channel has a cross-spectral density (PSD) in the frequency band 0.5 Hz to 70 Hz on two channels and the other. Shows the difference from the most probable function defined in the order of the q exponential function in the channel.
  • MF can be calculated by the formula (30) and shows the median value of the cumulative sum of MEG PSDs in each frequency region.
  • SE can be calculated by Eq. (31) and indicates the measure at which a certain spectral density is observed for each channel.
  • SSE can be calculated by Eq. (32) and indicates a measure in which a certain spectral density is observed in the entire set of MFs, which is a frequency that is 50% of the cumulative sum of each frequency band.
  • Metrics TSE, QSE, MF, SE, SSE can be calculated using spectral density PSD.
  • FIG. 13 is a schematic diagram showing an example of the magnetic field phase-locked distribution in the head of a dementia patient and a healthy person using a cross spectrum.
  • an Alzheimer's disease patient (AD) is given as an example of a dementia patient.
  • AD Alzheimer's disease patient
  • the difference in the cross-spectral distribution means the difference in the activity distribution of the cranial nerves. It can be seen that in the case of healthy subjects, the cross-spectrum values are higher in each part of the brain than in dementia patients.
  • a high cross-spectral value means that the activity of the cranial nerves is strongly activated.
  • FIG. 14 is a schematic diagram showing an example of the entropy distribution of subjects by age.
  • the horizontal axis represents age and the vertical axis represents metric SE.
  • Metrics SE can be calculated using the power spectrum PSD. It can be seen that the entropy (black dot) indicating a healthy person and the entropy (X mark) indicating a dementia patient show a stochastic difference.
  • This TSE uses equation (34), which expresses that it is composed of a set of entropy elements (escort distribution).
  • the entropy distribution consisting of the aggregate elements of entropy is unknown.
  • the equation (35) is used, which expresses that statistical analysis based on a probability distribution using the q order is possible by assuming a specific assumption.
  • Equation (35) determines the value of q that determines the probability distribution to which the entropy set belongs.
  • Equation (35) represents Jackson's q derivative operator in Tsallis statistics.
  • the test is possible by obtaining the abnormal probability using the Mahalanobis distance with respect to the probability distribution using the q order as an index.
  • the Mahalanobis distance means the information quantity distance when the q-Gaussian distribution shows a 0 mean multivariate distribution and shows an unbiased property.
  • the q-Gaussian distribution is a probability distribution obtained by maximizing the Tsallis entropy under constant dispersion conditions.
  • the brain state estimation unit 60 estimates the brain state of the subject based on the degree of deviation of the measure when the entropy is converted assuming that the entropy calculated by the entropy calculation unit 59 has a certain probability measure. can do.
  • the probability distributions are, for example, equilibrium probability distribution, normal distribution, discrete uniform distribution, simultaneous distribution, negative binomial distribution, geometric distribution, binomial distribution, Poisson distribution, Cauchy distribution, Including gamma distribution and so on.
  • the brain state estimation unit 60 estimates the brain state of the subject based on the degree of deviation of the entropy probability set when the calculated entropy probability distribution is assumed to be the equilibrium probability distribution.
  • the brain state estimation unit 60 can estimate the brain state by expressing the degree of divergence as the Mahalanobis distance of the brain disease group with respect to the healthy subject group.
  • the brain state estimation unit 60 can also discriminate the stimulation of the five senses, which is limited to healthy subjects and patients with brain diseases.
  • the brain state estimation unit 60 can be composed of a trained model generated by machine learning.
  • FIG. 15 is a schematic diagram showing an example of the relationship between the Mahalanobis distance and the anomaly probability.
  • the Mahalanobis distance is known to follow a chi-square distribution with m degrees of freedom.
  • ⁇ ⁇ (x- ⁇ ) T ⁇ -1 (x- ⁇ ) with the 95% point ⁇ of the chi-square distribution as a threshold value x can be determined as an abnormal value.
  • x is data sampled from the multivariate normal distribution N ( ⁇ , ⁇ ).
  • N ⁇ , ⁇
  • FIG. 16 is a schematic diagram showing the statistical properties (unbiasedness) of entropy sampled from each channel.
  • AD Alzheimer's disease
  • NC healthy subjects
  • UMVUs unbiased estimators
  • FIG. 17 is a schematic diagram showing an example of the results of outlier tests using entropy sampled from each channel.
  • the solid line shows the distance (measurement) between the two probability measures of NC with respect to NC
  • the broken line shows the distance (measurement) between the two probability measures of AD with respect to NC.
  • the degree of deviation of the probability set of entropy calculated individually for each MEG channel is expressed by the Mahalanobis distance. As can be seen from the figure, it can be verified that the distinction between the AD group and the NC group is statistically valid.
  • This embodiment includes dementia (including AD, DLB, frontal temporal lobe degeneration (FTLD), normal pressure hydrocephalus (NPH), etc.), brain tumor, psychiatric disorder (also referred to as mental illness, schizophrenia, epilepsy). , Mood disorder, addiction disorder, higher dysfunction, etc.), Parkinson's disease, Asperger's syndrome, attention defect / hyperactivity disorder (ADHD), sleep disorder, childhood disease, ischemic brain disorder, mood disorder (depression, etc.) It can be used for screening, differentiation, diagnosis, and classification of brain diseases including).
  • the display panel 61 can be composed of a liquid crystal panel, an organic EL (Electro Luminescence) display, or the like. In addition to or in place of the display panel 61, an external display device may be provided.
  • an organic EL Electro Luminescence
  • the operation unit 62 can be configured by, for example, a keyboard, a mouse, or the like.
  • the operation unit 62 may be composed of a touch panel incorporated in the display panel 61. The user can perform a predetermined operation on the display panel 61.
  • the display processing unit 63 performs a process of displaying information such as an estimation result output by the brain state estimation unit 60 on the display panel 61. As illustrated in FIG. 13, the display processing unit 63 can display the cross-spectral distribution between two channels among the plurality of channels corresponding to the plurality of positions in the brain.
  • the display processing unit 63 can display the cross-spectral distribution between two channels among the plurality of channels corresponding to the plurality of positions in the brain.
  • FIG. 18 is a schematic diagram showing a first example of the brain state estimation result.
  • the patient name, diagnosis result, and reliability are displayed on the display panel 61.
  • the diagnosis results indicate any of, for example, normal (healthy), mild cognitive impairment (MCI), Alzheimer's disease (AD), and other brain disorders.
  • the reliability indicates the certainty of the diagnosis result calculated based on the above-mentioned abnormality probability or the like.
  • An image in which a cross-spectral distribution is superimposed on a brain volume model can be displayed on the display panel. The displayed image can be rotated in a desired direction by the "rotate” icon, and the image can be enlarged or reduced by the "enlarge / reduce” icon.
  • the display processing unit 63 can superimpose the cross spectrum distribution on the brain volume model and display it.
  • the brain volume model can be constructed from MRI images, MI images, CT images and the like.
  • FIG. 19 is a schematic diagram showing a second example of the brain state estimation result.
  • a typical or typical cross-spectral distribution of a brain disease patient is displayed, and the cross-spectral distribution of the subject himself / herself is displayed in a contrastable manner.
  • the cross-spectral distributions of NC (healthy subjects), MCI (mild cognitive impairment), AD (Alzheimer's disease), and other brain diseases are displayed as patients with brain diseases.
  • the diagnosis result of the subject is MCI
  • the degree of dissociation (probability of abnormality probability) of the brain state of the brain disease group (MCI, AD, etc. in the example of the figure) with respect to the healthy subject group is It is displayed in contrast.
  • the display processing unit 63 can display the cross-spectral distributions of the healthy subject and the patient having a brain disease. In addition, the display processing unit 63 can display the degree of dissociation of the brain state of the brain disease group with respect to the healthy subject group.
  • FIG. 20 is a flowchart showing an example of the procedure of the brain state estimation process by the brain state estimation device 50.
  • the control unit 51 acquires the brain activity data of the subject (MEG data of each MEG channel) (S11).
  • MEG data data such as ocular potential (EOG), electroencephalogram (EEG), electrocardiogram (ECG), and myoelectric potential (EMG) can also be acquired.
  • EEG ocular potential
  • EEG electroencephalogram
  • ECG electrocardiogram
  • EMG myoelectric potential
  • the control unit 51 performs preprocessing including independent component analysis (ICA) on the acquired MEG data (S12).
  • Preprocessing other than ICA includes, for example, bias correction that subtracts the average value of MEG data from the acquired MEG data, a bandpass filter of a predetermined frequency band, a notch filter of a predetermined frequency, and an addition average for each predetermined period. At least one of can be used.
  • the control unit 51 performs SSP (signal space projection method) processing (S13), estimates the signal original position (S14), and calculates the signal distribution (S15).
  • the signal distribution is a magnetic field distribution.
  • the control unit 51 uses the MEG data obtained as a result of the processing in steps S14 or S15 to perform spectral analysis between each channel and between channels (S16). Here, the power spectrum and the cross spectrum are calculated.
  • the control unit 51 calculates the entropy (S17), and uses the calculated entropy to estimate the brain state of the subject based on the degree of deviation of the entropy probability set when the probability distribution of the entropy is assumed to be the equilibrium probability distribution. (S18), the estimation result is output (S19), and the process is terminated.
  • the brain state estimation device 50 can also be realized by using a computer equipped with a CPU (processor), RAM, and the like.
  • a recording medium on which a computer program for which a processing procedure is defined as shown in FIG. 20 is recorded is read by a recording medium reader, the read computer program is loaded into a RAM provided in the computer, and the computer program is loaded by a CPU (processor). By executing this, the brain state estimation device 50 can be realized on a computer.
  • MEG data has been described as an example of the subject's brain activity data, but the brain activity data is not limited to the MEG data, and the EEG data is also the present embodiment. Can be applied.
  • the brain state estimation device 50 can acquire MEG data from a biomagnetic measurement device (MEG measurement device) via wireless communication or wired communication.
  • MEG measuring device arranges a large number of MEG sensors (detection coils) on the scalp of the subject and outputs time-series data recorded by the MEG sensor. Channels are associated with each of the many MEG sensors.
  • the brain state estimation device of the present embodiment includes an acquisition unit that acquires brain data of a subject that fluctuates in a space including a time axis, and a generation unit that converts the acquired brain data to generate information related to the brain. Estimate to estimate the brain state of the subject based on the degree of divergence between the calculation unit that calculates the entropy based on the generated information and the measured entropy when the calculated entropy is assumed to have a certain probability measure. It has a part.
  • the acquisition unit obtains brain data including brain activity data, or brain data including brain activity data, myoelectric potential data, electrocardiogram data, electrooculogram data, and reference noise data.
  • the generation unit generates information related to the brain including frequency band information of rhythm or fractal information by using frequency conversion or conversion by nonlinear time series analysis, and the calculation unit generates the frequency of the generated rhythm.
  • the entropy that characterizes the amount of synchronous phase difference in the brain based on the information including the band information, or the entropy that characterizes the disorder of the time series signal and the degree of disruption of the phase structure is calculated based on the information including the generated fractal information.
  • the estimation unit estimates the brain state of the subject based on the degree of deviation of the entropy probability set when the calculated entropy probability distribution is assumed to be the equilibrium probability distribution.
  • the acquisition unit acquires time-series brain activity data for each of a plurality of channels corresponding to a plurality of positions in the brain, and the generation unit acquires each of the plurality of channels.
  • the frequency band information of the rhythm related to the brain activity is generated based on the brain activity data of the above and the brain activity data between two channels out of the plurality of channels.
  • the calculation unit calculates the entropy for each channel based on the measure at which a certain spectral density is observed.
  • the calculation unit calculates entropy based on a measure in which a certain spectral density is observed simultaneously in two channels.
  • the estimation unit estimates the brain state by expressing the degree of divergence as the Mahalanobis distance of the brain disease group with respect to the healthy subject group.
  • the brain state estimation device of the present embodiment includes a bias correction that subtracts the average value of the brain activity data from the brain activity data acquired by the acquisition unit, a bandpass filter of a predetermined frequency band, a notch filter of a predetermined frequency, and the like.
  • a first noise removing unit for removing noise of the brain activity data is provided by using at least one of the added averages for each predetermined period.
  • the brain state estimation device of the present embodiment is a second noise removing unit that removes noise from the brain activity data acquired by the acquisition unit by using independent component analysis and a signal space projection method. To be equipped.
  • the brain state estimation device of the present embodiment includes a display processing unit that displays a cross spectrum distribution between two channels among a plurality of channels corresponding to a plurality of positions in the brain.
  • the display processing unit superimposes and displays the cross spectrum distribution on the brain volume model.
  • the display processing unit displays the cross-spectral distribution of each of a healthy person and a patient having a brain disease.
  • the display processing unit displays the degree of dissociation of the brain state of the brain disease group from the healthy subject group.
  • the computer program of the present embodiment acquires brain data of a subject that fluctuates in a space including a time axis, converts the acquired brain data to generate information related to the brain, and is based on the generated information.
  • the entropy is calculated, and the process of estimating the brain state of the subject based on the degree of deviation of the measured entropy when the calculated entropy is assumed to have a certain probability measure is executed.
  • the brain state estimation method of the present embodiment acquires brain data of a subject that fluctuates in a space including a time axis, converts the acquired brain data to generate information related to the brain, and is based on the generated information.
  • the entropy is calculated, and the brain state of the subject is estimated based on the degree of deviation of the measured entropy when the calculated entropy is converted assuming that the calculated entropy has a certain probability measure.
  • the system of the present embodiment includes a function of classifying the brain activity data of a subject fluctuating in a space including a time axis by using a model in which learning data including a large number of solutions related to a brain function test is learned.
  • the system of the present embodiment uses the brain activity data of the subject, which fluctuates in the space including the time axis, by learning the solutions related to a large number of inputs and outputs including the activity area and the pattern of the brain. Includes the ability to estimate the spatial distribution of the signal source, including the time axis, as indicated by.
  • the system of the present embodiment uses a model that learns the spatial distribution including the time axis of the signal source estimated from the brain activity data of the subject that fluctuates in the space including the time axis and the case, and determines the brain state of the subject. Includes the ability to estimate.
  • the model includes a time series neural network.
  • the brain activity data includes brain activity data for a stimulus including at least one of visual, auditory, tactile, taste, and olfactory senses.
  • the method of the present embodiment includes classifying the brain activity data of a subject that fluctuates in space including the time axis by using a time-series model in which learning data including a large number of solutions related to a test of brain function is learned. ..
  • the brain activity data is obtained by using a model in which solutions related to a large number of inputs and outputs including brain activity regions and patterns are learned from the brain activity data of a subject that fluctuates in a space including a time axis. Includes estimating the spatial distribution including the time axis of the signal source indicated by.
  • the brain state of the subject is determined by using a model in which the spatial distribution including the time axis of the signal source estimated from the brain activity data of the subject fluctuating in the space including the time axis and the case are learned. Including estimating.
  • the computer program of the present embodiment has an instruction for operating as the above-mentioned system.
  • Brain state estimation device 51
  • Control unit 52
  • Data acquisition unit 53
  • Storage unit 54
  • Preprocessing unit 55
  • SSP processing unit 56
  • Signal source position estimation unit 57
  • Signal distribution calculation unit 58
  • Spectrum analysis unit 59
  • Entropy calculation unit 60
  • Brain condition estimation unit 61
  • Display panel 62
  • Operation unit 63 Display processing unit

Abstract

脳活動データから脳機能を検査するシステム、方法、脳状態推定装置、コンピュータプログラム及び脳状態推定方法を提供する。　時間軸を含む空間で変動する被験者の脳活動データを、脳機能の検査に係る多数の解を含む学習データを学習した時系列モデルを用いて分類する機能を含むシステムを提供する。このシステムは、脳活動データから、脳の活動領域およびパターンを含む多数の入出力に係る解を学習したモデルを用いて、前記活動データにより示される信号源の時間軸を含む空間分布を推定する機能を含んでもよく、推定された信号源の時間軸を含む空間分布と症例とを学習したモデルを用いて、被験者の脳状態を推定する機能を含んでもよい。

Description

脳状態推定装置、コンピュータプログラム、脳状態推定方法、脳機能の検査システムおよび方法
　本発明は、脳状態推定装置、コンピュータプログラム、脳状態推定方法、脳機能の検査システムおよび方法に関するものである。
　刺激などに対する脳の状態を、外部で得られる信号を介して捉えて、認知症の可能性などを含めた被験者（受験者、患者）の脳機能を検査することが試みられている。
　脳機能検査には一般に、被験者に刺激を提示し、それに対する脳の反応を測定する方法が用いられる。検査項目に応じ適切な刺激、測定、解析を、適切なパラメータを選択して実行し、さらに、得られたデータに対して適切な解を求めることが要望されている。
　本発明の一態様は、時間軸を含む空間で変動する被験者の脳活動データを、脳機能の検査に係る多数の解を含む学習データを学習した時系列モデルを用いて分類する機能を含むシステムである。時系列モデルの一例は、リカレントニューラルネットワーク（ＲＮＮ）を含む、時系列ニューラルネットワークである。他の態様の１つは、時間軸を含む空間で変動する被験者の脳活動データから、脳の活動領域およびパターンを含む多数の入出力に係る解を学習したモデルを用いて、活動データにより示される信号源の時間軸を含む空間分布を推定する機能を含むシステムである。他の態様の１つは、時間軸を含む空間で変動する被験者の脳活動データから推定された信号源の時間軸を含む空間分布と症例とを学習したモデルを用いて、被験者の脳状態を推定する機能を含むシステムである。これらのモデルは時系列ニューラルネットワークを含んでもよい。脳活動データは、安静時のデータであってもよく、視覚、聴覚、触覚、味覚、嗅覚の少なくともいずれかを含む刺激に対する脳活動データを含んでいてもよい。本発明は、コンピュータを上記のシステムとして稼働するための命令を有するプログラムを含む。
　本発明のさらに異なる他の態様の１つは、時間軸を含む空間で変動する被験者の脳活動データを、脳機能の検査に係る多数の解を含む学習データを学習した時系列モデルを用いて分類することを含む方法である。この方法は、時間軸を含む空間で変動する被験者の脳活動データから、脳の活動領域およびパターンを含む多数の入出力に係る解を学習したモデルを用いて、前記活動データにより示される信号源の時間軸を含む空間分布を推定することを含んでもよく、時間軸を含む空間で変動する被験者の脳活動データから推定された信号源の時間軸を含む空間分布と症例とを学習したモデルを用いて、被験者の脳状態を推定することを含んでもよい。
ＭＥＧの計測方法と賦活領域の位置推定の技術を示す図。 鑑別実験結果を示す図。 ＥＥＧ／ＭＥＧの周波数帯の特徴を参考に示す図。 本実施の形態の脳状態推定装置の構成の一例を示すブロック図である。 エポックデータの構成を示す模式図である。 加算平均処理の一例を示す模式図である。 ＭＥＧセンサマップ中心位置に対して頭部の回転に伴う脳表面活動分布の局所性を示す模式図である。 ＭＥＧセンサと電流双極子との間で構成される磁場分布の算出方法を示す模式図である。 本実施の形態により算出した磁場分布（Source Space Data）の結果の一例を示す模式図である。 周期性のずれが異なる二標本の時間領域での平均結果の一例を示す模式図である。 時間領域Ｔを用いて計算された自己共分散関数のＰＳＤの一例を示す模式図である。 認知症患者と健常者でのＰＳＤの違いの一例を示す模式図である。 クロススペクトルを用いた認知症患者と健常者の頭部内磁場位相同期分布の一例を示す模式図である。 年齢別の被験者のエントロピー分布の一例を示す模式図である。 マハラビノス距離と異常確率の関係の一例を示す模式図である。 各チャネルからサンプリングされたエントロピーの統計的性質（不偏性）を示す模式図である。 各チャネルからサンプリングされたエントロピーを用いた異常値検定の結果の一例を示す模式図である。 脳状態の表示画面の第１例を示す模式図である。 脳状態の表示画面の第２例を示す模式図である。 脳状態推定装置による脳状態推定処理の手順の一例を示すフローチャートである。
　一例として、脳磁図を用いたニューラルネットワークによる脳活動の分類を説明する。脳磁図（Magnetoencephalography、ＭＥＧ）は、脳磁場計測装置（脳磁計と略称される）を用いて記録される脳内の磁場活動である。幅広く普及している脳電位（Electroencephalography、ＥＥＧ）と対比して、本質的には同一現象を異なった方法で検索するものである。大脳皮質錐体細胞の尖端樹状突起のある部分が興奮して脱分極が生じると細胞外及び細胞内に電流が流れる。この細胞外電流を記録したものがＥＥＧであるのに対して、細胞内電流を取り巻くように生じる磁場を記録したものがＭＥＧである。ＭＥＧはＥＥＧと同様に脳から出る反応を記録するわけであり、脳に直接刺激あるいは負荷をかける必要が無いため、極めて安全な検査法であると言える。
　ＥＥＧと比べた場合のＭＥＧの最も大きな長所は高い空間分解能である。脳と頭皮の間には脳脊髄液、頭蓋骨、皮膚という導電率が大きく異なる３つの層がある。従って、脳で発生した電場はそれらによって大きな影響を受けるため、頭皮上に置いた脳波電極から脳の活動部位を正確に推測することは、特殊な推定法（例えば双極子追跡法）を用いない限り困難である。しかしながら、ＭＥＧの場合、磁場は導電率の影響を受けないため、記録条件が良好ならばｍｍ単位の精度で活動部位を推測することができる。
　ＭＥＧ、ＥＥＧと同様に脳の賦活状態とその領域を計測する手法にｆｕｎｃｔｉｏｎａｌ　ＭＲＩ（ｆＭＲＩ、機能画像）があるが、ｆＭＲＩと比べた場合に、ＭＥＧの長所としては、局所脳血流の変化ではなく細胞の電気的活動そのものを直接的に記録すること、及びＥＥＧと同様にミリ秒単位の高い時間分解能を有することが上げられる。
　本技術開発は、脳の活動状態が時間を含む空間で変動するデータとして得られる場合に、時系列モデルを用いて分類、推定および解析を行うことを提案する。ＭＥＧは、脳の活動データを得る方法として現在最適であるが、ＥＥＧにより得られたデータであってもよい。例えば、ＭＥＧの解析にニューラルネットワークを用い、脳の賦活情報を元に各種脳疾患の状態や健常人の脳活動の状態を高精度に鑑別するものであり、また信号源推定のような解析前処理へもニューラルネットワークを用いた方法を作成するものである。
　図１に示すように、ソースレベル（Source Level）データを時系列ＮＮ（ニューラルネットワーク）で分類するプロセス１００を設けてもよい。時系列ＮＮには、リカレントニューラルネットワーク等が含まれる。また、解のモデルを複数大量に用意しておき、入出力を再現するネットワークを学習する新規のプロセス２００を含む信号源推定方法を提供してもよい。また、信号源推定データを、ＮＮを用いて鑑別するプロセス３００を設けてもよい。鑑別する対象は、脳機能を含み、脳の異常に関する症状の有無を鑑別してもよい。例えば、認知症においては、ＮＣ／ＭＣＩ／ＡＤの分類等を行ってもよい。脳活動データは、安静時のデータであってもよく、視覚、聴覚、触覚、味覚、嗅覚の少なくともいずれかを含む刺激に対する脳活動データであってもよい。
　図２に、健常者の安静時データを、ＲｅｓＮｅｔモデルを用いて、自発閉眼と、自発開眼とにクラス分けした実験例を示している。このように、開眼閉眼に関する脳活動データを学習したモデルを用いることにより、脳活動データから眼の開閉に関するデータを抽出して解析できる。
　図３に、参考として、目の開閉に関するＥＥＧの周波数帯毎の電波強度を示し、また、健常者と認知症患者の周波数帯毎の電波強度に相違があることを示す。したがって、時間に依存して変動する脳活動データを取得し、適切な周波数帯に対する知識と、発生源に対する知識とを学習した人工知能を用いることにより、脳の機能の健全性または脳の機能の経時変化を、患者毎に的確に判断する方法およびシステムを提供できる。本発明における、時間を含む空間は、時間領域の関数をフーリエ変換などにより周波数空間に変換した空間も含む概念であり、静的なイメージあるいは信号ではなく、動的なイメージあるいは信号により得られる脳活動データを取り扱うシステムおよび方法に関するものである。
　脳内の磁場変化や電位変化を計算することで、電気生理学的活動領域で構成されたネットワーク（「機能的接続性」ともいう）を計算することができる。機能的接続性とは、神経細胞集団が構成する脳内のネットワークによって生じる加算された電気現象（ＥＰＳＰ：興奮性シナプス後電位、ＩＰＳＰ：抑制性シナプス後電位）の変化で表現される。解剖学的な神経投射経路の領域ごとに、入出力する磁場ベクトル強度または電位インピーダンスベクトルを計算した結果が機能的接続性になる。機能的接続性は、神経細胞集団に投射された電気的活動の経時的な情報を元に、統計的な関係性を表現する脳の接続性の分類の一種である。機能的接続性は、脳機能という観点から脳全体を捉えることで、脳内の異なる領域間の関係性をネットワークとして表す。機能的接続性を計算するには、まず生体磁気計測装置（ＭＥＧ計測装置）のＭＥＧセンサ（磁束検出コイル）から計測された磁場を用いて、脳表面の活動分布を推定するソースレベル解析（特に最小ノルム推定法(ＭＮＥ)）が使用されている。ＭＥＧ解析において被験者が安静した状態で観測される脳磁場分布から計算された機能的接続性を、デフォルトモードネットワーク（ＤＭＣ)という。ＭＮＥを用いたＤＭＣは、異なる脳領域で同期的に分散した神経活動分布を意味する。以下、脳状態の推定方法（健常者との差異を自動分類する方法）について説明する。
　図４は本実施の形態の脳状態推定装置５０の構成の一例を示すブロック図である。脳状態推定装置５０は、装置全体を制御する制御部５１、データ取得部５２、記憶部５３、前処理部５４、ＳＳＰ処理部５５、信号源位置推定部５６、信号分布算出部５７、スペクトル解析部５８、エントロピー算出部５９、脳状態推定部６０、表示パネル６１、操作部６２、及び表示処理部６３を備える。
　制御部５１は、ＣＰＵ（Central Processing Unit）、ＲＯＭ（Read Only Memory）及びＲＡＭ（Random Access Memory）などで構成することができる。記憶部５３は、半導体メモリ又はハードディスク等で構成することができる。以下では、時間軸を含む空間で変動する被験者の脳データとして、生体磁気計測装置のＭＥＧセンサ（検出コイル）によって計測される磁場データ（以下、「ＭＥＧデータ」とも称する）を例として説明する。脳データは、脳活動データ、眼電位データ（ＥＯＧ）、心電図データ（ＥＣＧ）、筋電位データ（ＥＭＧ）、基準ノイズデータ（reference noise）等を意味し、脳データは、脳活動データを含む場合、あるいは、脳活動データに加えて、ＥＯＧ、ＥＣＧ、ＥＭＧ、及びreference noiseのうちの少なくとも１つのデータを含む場合がある。例えば、脳状態を解析する場合には、脳活動データに加えて、ＥＯＧ、ＥＣＧ、ＥＭＧ、及びreference noiseのうちの少なくとも１つのデータを用いることができる。また、脳状態の差を用いて解析する場合には、脳活動データが必須となる。脳活動データは、磁場データに限定されるものではなく、ＥＥＧにより得られたデータでもよい。
　データ取得部５２は、時間軸を含む空間で変動する被験者の脳データを取得することができる。より具体的には、データ取得部５２は、脳活動データとしてのＭＥＧデータを取得することができる。生体磁気計測装置は、頭皮上に配置される多数（例えば、１６０か所など）のＭＥＧセンサを備え、各ＭＥＧセンサが各チャネルに対応しており、データ取得部５２は、複数のチャネルごとにＭＥＧデータを取得することができる。ＭＥＧデータは、例えば、被験者が安静した状態の臥位閉眼時のデータとすることができる。また、データ取得部５２は、ＭＥＧデータ以外のデータとして、眼電位（ＥＯＧ）、心電図（ＥＣＧ）、筋電位（ＥＭＧ）、基準ノイズデータなどのデータの全部又は一部を取得することができる。また、データ取得部５２は、ＭＥＧデータに加えて、あるいはＭＥＧデータに代えて、脳波（ＥＥＧ）を取得することができる。
　前処理部５４は、不良チャネルのデータ除去、ベースライン補正、ハンドパスフィルタ処理、ノッチフィルタ処理、加算平均処理、独立成分分析（ＩＣＡ：Independent Component Analysis）などを行う。計測された磁場データには、磁場由来のデータだけでなく、環境ノイズや頭部運動、生体ノイズ由来の信号が含まれている。前処理部５４は、磁場データ以外の信号を除去する機能を備える。
　不良チャネルのデータ除去は、計測できていないチャネルのＭＥＧデータを解析対象データから除去する。ＭＥＧデータの取得に失敗したチャネルの判定は、ＭＥＧデータ群に対する固有値分解によって求め、所定の閾値を示すチャネルのデータを除去する。除去したチャネルデータは、除去対象でないチャネルのデータで補間（内挿補間）することができる。
　ベースライン補正（「バイアス補正」とも称する）は、ＭＥＧデータの各チャネルに記録された時系列データの平均値を用いて、各チャネルに記録された時系列データから当該平均値を減算する。ＭＥＧデータは気温や気圧、個人差によって信号データの絶対値が変動するため、相対的なデータに変換するためである。
　ハンドパスフィルタ処理は、所要の周波数帯域をバンド幅とするバンドパスフィルタを用いて、ＭＥＧデータから不要な周波数帯域の信号を除去する。所要の周波数帯域は、例えば、０～７０Ｈｚとすることができる。この周波数帯域には、α波（８～１３Ｈｚ）、β波（１４～３０Ｈｚ）、θ波（４～７Ｈｚ）、δ波（０．５～３Ｈｚ）、γ波（２６～７０Ｈｚ）などの神経細胞集団から構成された脳の電気生理学的な電気信号が含まれる。
　ノッチフィルタ処理は、ノッチフィルタを用いて、生体磁気計測装置を稼動する際に使用する他の装置等の電源から混入する電源由来周期ノイズ（５０Ｈｚ又は６０Ｈｚ）を除去する。
　加算平均処理は、ＭＥＧデータに含まれるノイズ成分を除去するため、記録された時系列データを、例えば、５秒間の区間に分割し(trial)、分割した各データの加算平均を行うことでノイズ成分を減衰させる（エポックデータ）。
　図５はエポックデータの構成を示す模式図である。図５に示すように、複数のチャネルそれぞれのＭＥＧデータを複数の区間に分割する（trial１、trial２、…、trialｎ）。各区間の時間は、例えば、５秒間とすることができるが、これに限定されるものではない。また、区間の数ｎは適宜決定することができる。
　図６は加算平均処理の一例を示す模式図である。加算平均処理は、１つのtrial毎に１～ｎの各trialの信号を加算する。各trialに共通して観測される信号は加算されることにより強調される。一方、各trialでランダムに観測される波形（例えば、ノイズの波形）は、加算されることによって相殺され、ノイズの波形は低減する。
　独立成分分析は、単一の信号を、それが生じている独立な信号源ごとに分離する手法である。分離された信号に、ＭＥＧデータ、環境ノイズ、ＥＭＧデータが含まれる。独立成分分析は、室内環境に由来する環境ノイズを計測しているリファレンスチャネルと、眼球運動や筋肉の活動を計測しているＥＭＧチャネルの二つの情報を用いて、頭部に位置している各ＭＥＧセンサで計測したデータに混入されているノイズ成分（例えば、ＥＭＧ、リファレンスに同期しているノイズ）を除去する。具体的には、単一の被験者に対して計測された磁場データから、ランダム信号や頭部の運動に伴うアーティファクト、生体ノイズ由来の信号を除去し、複数の信号成分を求める。相関係数又は閾値に基づいて、アーティファクトとして定義するＥＭＧ－ＥＭＧ、ＥＭＧ－ＭＥＧの組み合わせのＩＣＡ要素を各チャネルから除去する。
　ＳＳＰ処理部５５は、ＳＳＰ（信号空間投影法： Signal Space Projection）手法を用いてＭＥＧデータから、信号源のノイズ影響を低減する。以下、ＳＳＰの方法について説明する。
　各ＭＥＧセンサで計測された信号ｙ（ｔ）は、元信号ｙｓ（ｔ）とノイズεの影響を含んでいるものと仮定すると、式（１）で表すことができる。
Figure JPOXMLDOC01-appb-M000001
　この時、時刻ｔにおける元信号ｙｓ（ｔ）は、様々な信号源の（回転、位置が不定の信号源）とその信号源に対して生成されるリードフィールドｌによって構成される。ある信号源Ｓｑは、ＭＥＧのチャネルに対して一意の分布を示し、式（２）で表される。ここで、ｑはチャネルインデックス、Ｑはチャネル数、ｌｑはチャネルｑの信号状態を示す。信号源Ｓｑは未知である。元信号ｙｓ（ｔ）とノイズεの影響にすべてに影響を及ぼすノイズＰは、式（３）で表すことができる。ノイズの影響Ｐを最小化した場合、式（１）から式（４）を定義できる。
　ノイズの影響Pの最小化には、式（５）で示す、Ｕｑの正規直交基底を考え、式（６）のように信号源Ｓｑ＝Ｃｑ、リードフィールドｌｑ＝Ｕｑと置き換えることで実現できる。データサンプル数が有限の場合、Ｕｑの正規直交基底からなる基底ベクトルｕは、データ共分散行列のＱ個の要素からなる最大固有値に対応した固有ベクトルの最尤推定解となる。式（６）を求めるために、データの共分散行列から特異値分解を用いることで、元信号ｙｓ（ｔ）を求めることができる。この元信号ｙｓ（ｔ）をセンサ空間データ（Sensor Space Data）とする。
　磁場データに、理論電圧値から構成した磁場分布に独立な乱数データと独立な環境で計測されたＥＭＧ＋ＥＣＧデータを混入し、上述のＩＣＡ処理及びＳＳＰ処理を行った。計測値と理論値の一致度合いを評価する指標としてＧＯＦ（Goodness of Fit）の計算結果は９７％と高かった。ＧＯＦは、ＧＯＦ＝｛１-Σ（計測値-理論値）2 ／Σ（計測値2 ）｝×１００［%］で求めることができる。ここで、Σは、各センサについて総和していることを意味する。ＧＯＦは、１００％のとき両者は完全に一致しており、１００％より小さな値になるほど一致度が小さくなることを意味する。なお、ＩＣＡ処理のみの場合、ＧＯＦは６０％と低かった。
　上述のように、本実施形態では、各チャネルの時系列データの無相関な信号や、チャネル間の無相関な信号を除去した上で、ＭＥＧデータを再構成することができる。出願人は、健常者であれば、脳内の機能的接続性が関連性を有し、同期性を有するのに対し、認知症などの脳疾患を有する患者では、機能的接続性の同期性が徐々に失われていくという点に着目した。本実施の形態によれば、チャネル間の無相関な信号を除去することにより、チャネル間の信号の相関成分を強調し、再構成されたＭＥＧデータを用いることにより、健常者と脳疾患患者とを分離（スクリーニング）することが可能になる。脳疾患の具体例としては、例えば、認知症（ＡＤ、ＤＬＢ、前頭側頭葉変性症（ＦＴＬＤ）、正常圧水頭症（ＮＰＨ）等を含む）、脳腫瘍、精神障害（精神疾患ともいう、統合失調症、てんかん、気分障害、依存障害、高次機能障害等を含む）、パーキンソン病、アスペルガー症候群、注意欠陥・多動性障害（ＡＤＨＤ）、睡眠障害、小児疾患、虚血性脳障害、気分障害（うつ病等を含む）等を含む。
　信号源位置推定部５６は、ＳＳＰ処理部５５で得られたセンサ空間データを用いて、脳内の信号源の位置を推定する。ここで、信号源は、多数のニューロンの同期活動を表す電流双極子である。磁場は距離の関数に依存する。ＭＥＧセンサの計測データはＭＥＧセンサと頭部間の距離の影響を強く受けている。
　図７はＭＥＧセンサマップ中心位置に対して頭部の回転に伴う脳表面活動分布の局所性を示す模式図である。図中、輪郭部分は、頭部中心位置及び範囲に従って回転拡縮された頭部範囲を示す。輪郭内のメッシュは、頭部中心位置及び範囲に従って回転拡縮されたＭＮＩテンプレートを示す。輪郭内の左右の部分は、頭部とＭＥＧセンサ位置の中で最も短い距離のエリアを示す。輪郭の周囲の点は、各ＭＥＧセンサの位置を表す。
　ＭＥＧセンサと頭部間の距離の影響を低減するには、頭部位置を考慮したセンサ位置を定義する必要がある。ＭＥＧデータに記録された鼻中央、右耳、左耳の空間位置から、頭部の中心位置を設定し、設定した中心位置に対応して各ＭＥＧセンサの位置が設定されるように、アフィン変換を用いて頭部モデル（ＭＮＩテンプレート）を回転拡縮させる。この頭部モデルの形状とセンサ距離をもとに、ＳＳＰ処理部５５で得られたセンサ空間データを使用して、媒質内を通過する磁場モデル式であるビオサバール式を使用して、脳内信号源位置を推定する（順問題によるダイポール推定）。
　信号源位置推定部５６は、推定した脳内信号源位置における電流成分から再びビオサバール式を用いて、頭部の外に位置するセンサまでの距離によって特異的に減衰、増強される磁場データに構成し、再構築されたセンサ空間データとする。
　信号分布算出部５７は、信号源位置推定部５６で再構築されたセンサ空間データを用いて、磁場分布（Source Space Data）を算出する。本実施の形態では、アダプティブビームフォーマー（ＳＡＭ）を用いて、磁場分布（Source Space Data）を算出する。
　図８はＭＥＧセンサと電流双極子との間で構成される磁場分布の算出方法を示す模式図である。ＭＥＧの磁場データは、電流成分を持つ双極子Ｊｏの脳内での位置と回転が与えられたとき(ｒ、ｒ-ｒ′)、外部で計測しているＭＥＧセンサに到達するまでの磁場の分布Ａ（ｒ）を図中の式であるビオサバール方程式を用いて算出することができる。ＳＡＭは、各ＭＥＧセンサで記録した磁場データから、ビオサバール方程式を用いて電流双極子Ｊｏの存在を推定した上で、数理計画法を用いて外部センサ（Sensor Space Data）、もしくは頭部表面の境界面（図中の円状領域）において計測される磁場分布（Source Space Data）を計算したものである。この手法を用いることで、計測された磁場データ内に含まれる磁場分布以外のシグナルを低減し、空間ＳＮＲを向上することができる。
　図９は本実施の形態により算出した磁場分布（Source Space Data）の結果の一例を示す模式図である。図に示すように、アダプティブビームフォーマー（ＳＡＭ）を用いて求められた磁場分布（Source Space Data）により、空間ＳＮＲが改善され、隣接した神経活動源を正確に評価することができる。
　次に、スペクトル解析について説明する。
　スペクトル解析部５８は、データ取得部５２を介して取得した脳データを変換して脳に関連する情報を生成することができる。例えば、スペクトル解析部５８は、周波数変換を用いて、律動の周波数帯域情報を含む脳に関連する情報を生成することができる。律動の周波数帯域情報は、脳機能情報が一定の周期構造を持っていることに注目した情報である。また、スペクトル解析部５８は、非線形時系列解析による変換を用いて、フラクタル情報を含む脳に関連する情報を生成することができる。フラクタル情報は、脳機能情報に含まれる乱雑性に注目した情報であり、脳機能情報の類似構造の破れ具合を意味する。具体的には、スペクトル解析部５８は、信号源位置推定部５６、又は信号分布算出部５７で得られた磁場データ（各ＭＥＧチャネルの時系列データ）を取得して、単一チャネルのスペクトル解析、及び２つのチャネル間のスペクトル解析を行うことができる。スペクトル解析の切り替えは、計算条件を含んだ設定ファイルを指定することで行うことができる。以下では、まず、単一チャネルのスペクトル解析について説明する。
　ウィンナー・フィンチンの定理（自己共分散関数のフーリエ変換はスペクトルであり、スペクトルの逆フーリエ変換は自己共分散関数である）より、式（７）で表す、自己共分散関数ρ(τ)のフーリエ変換を用いて、式（８）のように、スペクトルを解析する（相関法）。ここで、ｙ（ｔ）（ｔ＝０、１、…、Ｎ）は、信号列を示し、ｘ（ｔ）は、τインデックス分、離散的にサンプリングした信号列を示し、ｋはサンプリングインデックス数（ｋ＝０、１、…、Ｎ－１）を示す。
Figure JPOXMLDOC01-appb-M000002
　一般的に、スペクトルの推定は、有限のデータに対して行うために、データの端の影響が、推定されたスペクトルピークが存在する周波数から、他の周波数にスペクトルの漏れ(spectral leakage)として現れる問題がある。観測されたスペクトルに含まれるスペクトルの漏れにより、スペクトルは真の物理現象を完全に表現するものではないことが知られている。このため、窓関数の種類や、平滑化処理によってスペクトル推定の誤差を低減することができるが、推定されたスペクトルの誤差は大きいのが実情である。
　式（８）で示すような相関法を用いることで、任意に設定した時間領域で、ある時刻ｔからτだけ離れた各時点（ｔ＋τ）での関連性を定量的に評価することができる。すなわち、自己相関関数（Auto Correlation）からスペクトルを求めるBlackman-Tukey法は、スペクトルの誤差論に基づいて有限の時間領域での変動を求めるように設計されている。有限の時間領域におけるランダムシーケンスからなる集合要素は、各要素において個別の発生確率を有している。確率過程の構成要素である各時間領域の確率変数が独立に分布している場合には、異時点間の自己相関は全て０になる。ＭＥＧの時系列データに対して加算平均処理を行う場合、任意の時間領域を設定するため周期位相がずれた信号を同時に加算平均することになり、自己相関が小さく、ＤＴＦＴ（Discrete-time Fourier transform：離散時間フーリエ変換）を持たないことがある。
　図１０は周期性のずれが異なる二標本の時間領域での平均結果の一例を示す模式図である。図１０に示すように、複数の標本の属する確率過程を用いることによって統計量を推定できない場合は、標本として扱う時間領域を変更する必要がある。すべての統計量が任意の１つの標本過程ｘ（ｔ）から、一致性と不偏性を持って推定できる場合、すなわち、エルゴード性を満たす場合、統計量を示すことができる。従って、標本として扱う時間領域の変更によってエルゴード性を保証することを示さなければならない。
　本実施の形態では、平均μと自己共分散関数ρ（τ）（Autocovariance）についてのエルゴード性を保証する時間領域Ｔを自動決定する。時間領域での平均は、式（９）で表され、時系列データの平均エルゴード性は、式（１０）で表され、定常過程を示す分散の０収束条件は、式（１１）で表される。式（１１）において、Ｃx （ｔ１-ｔ２）は周波数領域ｔ１、ｔ２での自己共分散関数を示す。
Figure JPOXMLDOC01-appb-M000003
　特に、平均値に収束する性質と、分散が０に収束する条件の二つを満たす時間領域Ｔを探索し、その時間領域Ｔでの自己共分散関数ρ（τ）＝Ｐ（ω）を用いてパワースペクトル密度（ＰＳＤ）を推定することができる。自動決定した時間領域Ｔを用いたＰＳＤ計算は次のように行う。ある時間領域Ｔを持つ時系列データの矩形サンプリングを式（１２）とすると、この時のＰＳＤは、式（１３）で求めることができる。ここで、Ｅは式（１２）で規格化したτインデックス領域での時系列データｘτの、ｘτ（ｔ）ｘτ（ｔ＋τ）の時間平均を示す。τはサンプリング幅を示す。
Figure JPOXMLDOC01-appb-M000004
　図１１は時間領域Ｔを用いて計算された自己共分散関数のＰＳＤの一例を示す模式図である。図中、横軸は周波数を示し、縦軸はＰＳＤ（パワースペクトル密度）を示す。実線は理論値を示し、破線は時間領域ＴにおけるBlackman-Tukey法を用いた値を示す。矩形窓、バートレット窓、ブラックマン窓など、窓関数の種別によらず、最適な時間領域Ｔを求めることができる。
　図１２は認知症患者と健常者でのＰＳＤの違いの一例を示す模式図である。図中、横軸は周波数を示し、縦軸はＰＳＤ（パワースペクトル密度）を示す。実線は健常者（ＮＣ）を示し、破線はアルツハイマー病患者（ＡＤ）を示す。図に示すように、α帯、θ帯に特徴的な差が見られ、認知症の臨床段階の違いをスペクトルで定義することができる。逆に言えば、スペクトル密度を用いて認知症の臨床段階を分類、スクリーニングすることが可能となる。なお、図示していないが、認知症（ＡＤ、ＤＬＢ、前頭側頭葉変性症（ＦＴＬＤ）、正常圧水頭症（ＮＰＨ）、軽度認知障害（ＭＣＩ：Mild cognitive impairment）、アルツハイマー病による軽度認知障害（ＭＣＩdue to ＡＤ）、前駆期アルツハイマー病（prodromal ＡＤ）、アルツハイマー病の発症前段階／プレクリニカルＡＤ（preclinical ＡＤ）等を含む）、多発性硬化症脳腫瘍、精神障害（精神疾患ともいう、統合失調症、てんかん、気分障害、依存障害、高次機能障害、等を含む）、パーキンソン病、アスペルガー症候群、注意欠陥・多動性障害（ＡＤＨＤ）、睡眠障害、小児疾患、虚血性脳障害、気分障害（うつ病等を含む）等ついても分類、鑑別、診断、スクリーニングなどを期待できる。
　次に、２つのチャネル間のスペクトル解析について説明する。
　２つのチャネル間のスペクトル（「クロススペクトル」ともいう）解析には、脳活動を表現する７０Ｈｚ未満の周波数成分を持つ比較的短い時系列データに対して有効である、Multi-Taper Methodを使用したフーリエ変換を使用してクロススペクトルを解析する(直接法）。
　単一チャネルのスペクトル解析と同様に、クロススペクトルの推定は、有限のデータに対して行うために、データの端の影響が、推定されたスペクトルピークが存在する周波数から、他の周波数にスペクトルの漏れ(spectral leakage)として現れる問題がある。このスペクトルの漏れは、推定スペクトルのシグナル（平均）とノイズ（標準偏差）が同じ大きさを持ち、時間領域Ｔを変化させても推定スペクトルの誤差は改善されない場合がある。そこで、スペクトルの漏れを低減するために、テーパー(taper)関数を用いることができる。これをピリオドグラムの平滑化という。本実施の形態でのピリオドグラムの平滑化処理は、以下のとおりである。
　ＭＥＧの時系列データをｙ（ｔ）とし、テーパー関数をｈ（ｔ）とする。この時のフーリエ変換を行ったＹ′（ｋ）との関係は、式（１４）で表される。
Figure JPOXMLDOC01-appb-M000005
　フーリエ変換であるＹ′（ｋ）の逆フーリエ変換は、式（１５）であり、式（１５）を式（１４）に代入すると、式（１６）が得られる。式（１６）のＨ（ｋ-ｌ）は、テーパー関数ｈ（ｔ）のフーリエ変換式である、式（１７）で表すことができる。テーパーを用いたフーリエ変換の期待値{Ｓ′}は、時系列データｙ（ｔ）のフーリエ変換、及びテーパー関数ｈ（ｔ）のフーリエ変換との畳み込みで表される。フーリエ変換の期待値{Ｓ′}は、式（１８）で求めることができる。ここで、ｆN は解析対象である全周波数帯域を示し、Ｎ＝（０、１、…、Ｎ－１）［Ｈｚ］である。ｆ′はサンプリング周波数を示す。テーパー関数ｈ（ｔ）を用いたスペクトルデータは、式（１９）である。
　２つのチャネル間の組み合わせｊによるクロススペクトル推定は、式（１６）を拡張する式（２０）で表すことができる。２つのチャネル間の組み合わせｊ＝１、２の場合、クロススペクトル推定は、式（２０）により、式（２１）で表される。クロススペクトルを計算する２つのチャネル間のデータの位相が同期している場合、クロススペクトルは高い値を示す。
　テーパー関数を用いた修正ピリオドグラムの計算により、フーリエ窓関数（ハミング窓、矩形窓など）によらず、クロススペクトルの期待値は同じであり、クロススペクトル推定誤差の変動を極めて小さくすることができる。
　次に、認知症の特徴量としてのＰＳＩ（Phase Slope Index）について説明する。ＰＳＩは、ビームフォーマーを用いて脳磁場分布を定位させたメトリクス群の内、コヒーレンスを用いて計算したものであり、安定した統計量を示すことが報告されている。
　ｓk （ｔ）、ｓj （ｔ）を、ｋ点とｊ点での時系列磁場信号とすると、ｋ点とｊ点を結ぶある脳内のスペクトルを式（２２）、式（２３）で表すとする。
Figure JPOXMLDOC01-appb-M000006
　この時のクロススペクトルは、式（２４）で表すことができる。位相情報が第ｊ番目の位置から第ｋ番目の位置に伝わる時間をτとすると、クロススペクトルの同期位相状態は、信号の周波数ｆと時間差τに比例する。この情報をもとにPhase Slope（Ψ）は、式（２５）によって定義される。これは、クロススペクトルから計算できるが、コヒーレンスとしても記述できる。コヒーレンスは、クロススペクトルΦから、式（２６）として記述できる。コヒーレンスは、複素数であるクロススペクトルの実部を表現したもので、２つの信号のフーリエ周波数成分の相互相関係数と見ることができ、式（２６）で示す値は、０から１の間の値をとり、０に近い場合、両者の相関性は低く、１に近い場合、両者の相関性が高いことを意味する。ＰＳＩは、コヒーレンスを用いて、式（２７）によって求められる。式（２５）と式（２７）は同じとなる。ＰＳＩは、位相情報が第ｊ番目の位置から第ｋ番目の位置に伝わるラグを持った同期情報を意味する。
　前述のように、ＰＳＩを計算する際に必要となるコヒーレンスは、信号源近くにおいて、偽の磁場分布を形成することが知られている。このため、コヒーレンスを用いたＰＳＩの信頼性は低い。以下、本実施の形態では、コヒーレンスを使用せず、クロススペクトル、パワースペクトル計算を使用したエントロピーの算出について説明する。
　ＭＥＧデータを用いてスペクトル解析された結果、すなわち、スペクトルデータから脳活動に関連する各脳内領域間で協調した律動の周波数帯域（「同期成分、コヒーレンス」ともいう）、ある律動の周波数帯域でのパワー（「周波数成分」ともいう）の両方の情報を用いて、脳内の同期程度量（ある周波数帯域で見られるコヒーレンスが示す同期位相差量）と同期分布（コヒーレンスやクロススペクトルでの同期位相差量）を計算することで、脳疾患に関連した特徴量（脳疾患の臨床段階での差異、例えば、認知症患者と健常者の差異を決定付けるエントロピー）を算出することができる。
　エントロピー算出部５９は、スペクトル解析部５８が生成した情報に基づいてエントロピーを算出する。具体的には、２つのチャネル間に関連して、エントロピー算出部５９は、脳活動に関連する律動の周波数帯域情報に基づいて脳内の同期位相差量を特徴付けるエントロピーを算出することができる。周波数帯域情報は、律動の周波数帯域、及びある律動の周波数帯域でのパワー（「周波数成分」ともいう）を含む。同期位相差量は、ある周波数帯域で見られるコヒーレンスが示す同期程度量、コヒーレンスやクロススペクトルでの同期成分をいう。コヒーレンスは、グリッド間で同期している成分量を表すものであり、同じ位相、同じ周波数の信号成分をいう。また、各チャネルに関連して、エントロピー算出部５９は、脳活動に関連する律動の周波数帯域情報に基づいて脳内の同期量を特徴付けるエントロピーを算出することができる。また、各チャネルに関連して、脳活動に関連する時系列信号の構造に注目し、エントロピー算出部５９は、非線形時系列解析を用いて、時系列信号の乱雑性や位相構造の破れ具合を特徴付けるエントロピーを算出することができる。非線形時系列解析を用いて算出される、時系列信号の乱雑性や位相構造の破れ具合を特徴付けるエントロピーとしては、例えば、（１）信号の乱雑性指標、指定した組み合わせの時系列データの差分が、ある数値範囲でみられる測度を意味するapprocimate entropy（ＡｐＥｎ)、（２）自己相関を持つ信号の乱雑性指標を意味し、ＡｐＥｎの拡張系でデータサイズに依存しないSample Entropy、（３）コロモゴロフ過程（一様連続性）を仮定したときの時系列データの測度（乱雑性指標）を意味するLempel-Ziv complexity（ＬＺＣ）、（４）自己回帰係数を用いた時系列データの測度を意味するmultisacale entropy（ＭＳＥ）、（５）自己回帰係数を用いた時系列データの測度を意味するrefined ＭＳＥ（ＲＭＳＥ）、（６）ＬＺＣの拡張系であり、二つの時系列データの測度を比較するLopez Ruiz Mancini Calbet complexity（ＬＭＣ）などがある。
　以下では、単一チャネルに関するエントロピーを表すメトリクスＴＳＥ、ＭＦ、ＳＥ、ＳＳＥ、及び２つのチャネル間に関するエントロピーを表すメトリクスＱＳＥについて説明する。
　ＴＳＥは、式（２８）で計算することができ、各チャネルの０．５Ｈｚから７０Ｈｚの周波数帯域での、あるスペクトル密度（ＰＳＤ）が観測される測度と、ｑ指数関数のオーダーで定義された最尤関数との差分を示す。
Figure JPOXMLDOC01-appb-M000007
　ＱＳＥは、式（２９）で計算することができ、２つのチャネルで、０．５Ｈｚから７０Ｈｚの周波数帯域での、一方のチャネルでのあるクロススペクトル密度（ＰＳＤ）が観測される測度と他方のチャネルでのq指数関数のオーダーで定義された最尤関数との差分を示す。
　ＭＦは、式（３０）で計算することができ、ＭＥＧのＰＳＤの各周波数領域での累積和の中央値を示す。
　ＳＥは、式（３１）で計算することができ、各チャネルの、あるスペクトル密度が観測される測度を示す。
　ＳＳＥは、式（３２）で計算することができ、各周波数帯域の累積和の５０％となる周波数であるＭＦの全体集合の中で、あるスペクトル密度が観測される測度を示す。メトリクスＴＳＥ、ＱＳＥ、ＭＦ、ＳＥ、ＳＳＥは、スペクトル密度ＰＳＤを用いて計算することができる。
　図１３はクロススペクトルを用いた認知症患者と健常者の頭部内磁場位相同期分布の一例を示す模式図である。図では、認知症患者の例としてアルツハイマー病患者（ＡＤ）を挙げている。図から分かるように、健常者と認知症患者との間で、クロススペクトル分布に差異がある。別言すれば、クロススペクトル分布に基づいて、健常者と認知症患者の分類、スクリーニングを行うことができる。クロススペクトル分布の差異は、脳神経の活動分布の差異を意味する。認知症患者に比べて健常者の場合、脳内の各部においてクロススペクトルの値が高いことが分かる。クロススペクトルの値が高いということは、脳神経の活動が強く活性化していることを意味する。
　図１４は年齢別の被験者のエントロピー分布の一例を示す模式図である。図中、横軸は年齢を示し、縦軸はメトリクスＳＥを示す。メトリクスＳＥは、パワースペクトルＰＳＤを用いて計算することができる。健常者を示すエントロピー（黒点）と認知症患者を示すエントロピー（Ｘ印）とは、確率的な差異を示すことが分かる。
　次に、エントロピーを用いた検定について説明する。
　着目しているある物理量をＱとし、測定によって得られたその物理量の期待値をＱＥとする。しかし、それ以外の情報を取得することが困難である系では妥当な解析ができない。そこで、物理量Ｑについての事前確率分布を拘束条件として、系の状態に関する不確定性を最大化することで、その系の持つ性質を推定する。これを最大エントロピー原理といい、この要請に基づいて統計系の平衡確率分布が決定される。しかし、一般的な脳活動はマルチフラクタル系のような不均衡状態であり、従来の平衡した状態を解析する手法は妥当ではない。そこで、Boltzmann-Shannonエントロピーを一般化したTsallis統計を用いてマルチフラクタル系のように確率分布関数がベキ則的振る舞いをする場合に対応した統計系の平衡確率分布を計算し、この平衡確率分布を用いた指標計算を行うことができる。
　正規分布のような確率分布の場合、平均と分散が分かれば、分布に属する値の確率分布が分かる。しかし、エントロピーは限られたデータしか存在しないので、エントロピーの確率分布は分からない。そこで、以下のように、平衡分布（ｑガウス分布）とマハラビノス距離を用いることによって、エントロピーの値から健常者群と疾患群とを分類することができる。
　ＭＥＧチャネルのエントロピー分布に対して、Hotellingの不偏性検定を用いた多変量２標本問題の統計解析を可能にする。具体的には、以下のような手順を用いることができる。すなわち、
（１）脳内のエントロピーＴＳＥを表す式（３３）を使用する。
Figure JPOXMLDOC01-appb-M000008
（２）このＴＳＥはエントロピーの集合要素から成り立っていることを表す式（３４）を使用する（エスコート分布）。エントロピーの集合要素からなるエントロピー分布は未知である。
（３）未知のエントロピー分布に対して、特定の仮定を想定することにより、ｑオーダーを用いた確率分布に基づく統計解析が可能であることを表す式（３５）を使用する。式（３５）は、エントロピー集合の属する確率分布を決定づけるｑの値を決定する。式（３５）は、Tsallis統計におけるJacksonのｑ微分演算子を表す。
（４）ｑオーダーを用いた確率分布に対するマハラビノス距離を指標として異常確率を求めることにより、検定を可能とする。マハラビノス距離は、ｑガウス分布が０平均多変量分布を示し、不偏性を示す場合の情報量距離を意味する。ｑガウス分布は、分散一定の条件下でTsallisエントロピーを最大化して得られる確率分布である。
　上述のように、脳状態推定部６０は、エントロピー算出部５９が算出したエントロピーがある確率測度を持つと仮定した場合のエントロピーを変換した際の測度の乖離度に基づいて被験者の脳状態を推定することができる。ある確率測度を持つと仮定した場合の確率分布は、例えば、平衡確率分布、正規分布、離散一様分布、同時分布、負の二項分布、幾何分布、二項分布、ポワソン分布、コーシー分布、ガンマ分布などを含む。確率分布を平衡確率分布とした場合、脳状態推定部６０は、算出したエントロピーの確率分布を平衡確率分布と仮定した場合のエントロピーの確率集合の乖離度に基づいて被験者の脳状態を推定することができる。具体的には、一度ｑガウス分布に変換することで、マハラビノス距離という測度に対して属する確率分布を定義し、定義した確率分布間での乖離度を計算できる。これにより、脳状態推定部６０は、乖離度を健常者群に対する脳疾患群のマハラビノス距離で表現して脳状態を推定することができる。なお、脳状態推定部６０は、健常者や脳疾患患者に限らす、五感の刺激の判別を行うこともできる。脳状態推定部６０は、機械学習によって生成された学習済みモデルで構成することができる。
　図１５はマハラビノス距離と異常確率の関係の一例を示す模式図である。ｍ次元正規分布からデータをサンプリングしている場合、マハラノビス距離は、自由度ｍのカイ二乗分布に従うことが知られている。カイ二乗分布の９５％点αなどを閾値にして、α＜（ｘ-μ）T Σ-1 （ｘ-μ）のとき、ｘを異常値と判定することができる。ｘは、多変量正規分布Ｎ（μ、Σ）からサンプリングされたデータである。図１５の例では、例えば、確率が９５％に対応するマハラビノス距離より小さい場合、異常であると判定することができる。
　図１６は各チャネルからサンプリングされたエントロピーの統計的性質（不偏性）を示す模式図である。マハラビノス距離を用いて異常確率を求める場合、サンプリングされたエントロピーが統計的性質を持っている必要がある。図１６に示すように、認知症患者（ＡＤ：アルツハイマー病患者）と健常者（ＮＣ）の被験者間においても、全脳分布の各ＭＥＧチャネルからサンプリングされたエントロピーから計算された不偏推定量（ＵＭＶＵ）が一様の傾向を示すことが分かる。すなわち、エントロピーの分布は定常である。ここでは、エントロピーを表すメトリクスとしてＴＳＥを用いているが、他のメトリクスの場合も同様である。ＵＭＶＵ性質を持つ分布からサンプリングされたエントロピー分布の普遍性を示す統計的性質からカイ二乗分布として定義することができる。これにより、ＭＥＧチャネルのエントロピー分布はHotellingの不偏性検定を用いた多変量２標本問題の統計解析を可能にする。
　図１７は各チャネルからサンプリングされたエントロピーを用いた異常値検定の結果の一例を示す模式図である。実線は、ＮＣに対してＮＣの二つの確率測度との間の距離(計量)を示し、破線はＮＣに対してＡＤの二つの確率測度との間の距離(計量) を示す。各ＭＥＧチャネルで個別に計算したエントロピーの確率集合の乖離度をマハラノビス距離で表現したものである。図から分かるように、ＡＤ群とＮＣ群の区別が統計的に妥当であることが検証可能となる。本実施の形態は、認知症（ＡＤ、ＤＬＢ、前頭側頭葉変性症（ＦＴＬＤ）、正常圧水頭症（ＮＰＨ）等を含む）、脳腫瘍、精神障害（精神疾患ともいう、統合失調症、てんかん、気分障害、依存障害、高次機能障害等を含む）、パーキンソン病、アスペルガー症候群、注意欠陥・多動性障害（ＡＤＨＤ）、睡眠障害、小児疾患、虚血性脳障害、気分障害（うつ病等を含む）等を含む脳疾患のスクリーニング、鑑別、診断、分類に用いることができる。
　表示パネル６１は、液晶パネル又は有機ＥＬ（Electro Luminescence）ディスプレイ等で構成することができる。なお、表示パネル６１に加えて、あるいは代えて、外部の表示装置を設ける構成でもよい。
　操作部６２は、例えば、キーボード、マウスなどで構成することができる。操作部６２は、表示パネル６１に組み込まれたタッチパネルで構成してもよい。ユーザは表示パネル６１上で所定の操作を行うことができる。
　表示処理部６３は、脳状態推定部６０が出力する推定結果などの情報を表示パネル６１に表示する処理を行う。図１３に例示したように、表示処理部６３は、脳内の複数の位置に対応する複数のチャネルのうちの２つのチャネル間のクロススペクトル分布を表示することができる。以下、他の表示例について説明する。
　図１８は脳状態推定結果の第１例を示す模式図である。表示パネル６１には、患者名、診断結果、信頼度が表示される。診断結果は、例えば、正常（健常者である）、軽度認知障害（ＭＣＩ）、アルツハイマー病（ＡＤ）、その他の脳疾患などのいずれかを示す。信頼度は、前述の異常確率などに基づいて算出された、診断結果の確からしさを示す。表示パネルには、脳体積モデルにクロススペクトル分布を重畳した画像を表示することができる。表示された画像は、「回転」アイコンにより所望の向きに回転させることができ、「拡大・縮小」アイコンにより画像を拡大、縮小することができる。
　このように、表示処理部６３は、脳体積モデルにクロススペクトル分布を重畳して表示することができる。脳体積モデルは、ＭＲＩ画像、ＭＩ画像、ＣＴ画像などから構築することができる。
　図１９は脳状態推定結果の第２例を示す模式図である。表示パネル６１には、脳疾患患者の典型的又は代表的なクロススペクトル分布が表示されるとともに、被験者自身のクロススペクトル分布が対比可能に表示される。図の例では、脳疾患患者として、ＮＣ（健常者）、ＭＣＩ（軽度認知障害）、ＡＤ（アルツハイマー病）、その他の脳疾患それぞれのクロススペクトル分布が表示されている。また、図の例では、被験者の診断結果がＭＣＩであるとし、健常者群に対する脳疾患群（図の例では、ＭＣＩ、ＡＤ、その他）の脳状態の乖離度合い（異常確率の確からしさ）が対比して表示されている。
　このように、表示処理部６３は、健常者と脳疾患を有する患者それぞれのクロススペクトル分布を表示することができる。また、表示処理部６３は、健常者群に対する脳疾患群の脳状態の乖離度合いを表示することができる。
　次に、脳状態推定装置５０の処理について説明する。
　図２０は脳状態推定装置５０による脳状態推定処理の手順の一例を示すフローチャートである。以下では、便宜上、処理の主体を制御部５１として説明する。制御部５１は、被験者の脳活動データ（各ＭＥＧチャネルのＭＥＧデータ）を取得する（Ｓ１１）。なお、ＭＥＧデータ以外に、眼電位（ＥＯＧ）、脳波（ＥＥＧ）、心電図（ＥＣＧ）、筋電位（ＥＭＧ）などのデータを取得することもできる。
　制御部５１は、取得したＭＥＧデータに対して、独立成分分析（ＩＣＡ）を含む前処理を行う（Ｓ１２）。ＩＣＡ以外の前処理には、例えば、取得したＭＥＧデータに対してＭＥＧデータの平均値を減算するバイアス補正、所定の周波数帯域のバンドパスフィルタ、所定の周波数のノッチフィルタ、所定期間毎の加算平均の少なくとも１つを用いることができる。
　制御部５１は、ＳＳＰ（信号空間投影法）処理を行い（Ｓ１３）、信号原位置推定を行い（Ｓ１４）、信号分布を算出する（Ｓ１５）。ここでは、信号分布は、磁場分布である。制御部５１は、ステップＳ１４又はＳ１５の処理の結果得られた、ＭＥＧデータを用いて、各チャネル及びチャネル間のスペクトル解析を行う（Ｓ１６）。ここでは、パワースペクトル、クロススペクトルが算出される。
　制御部５１は、エントロピーを算出し（Ｓ１７）、算出したエントロピーを用いて、エントロピーの確率分布を平衡確率分布と仮定した場合のエントロピーの確率集合の乖離度に基づいて被験者の脳状態を推定し（Ｓ１８）、推定結果を出力し（Ｓ１９）、処理を終了する。
　脳状態推定装置５０は、ＣＰＵ（プロセッサ）、ＲＡＭなどを備えたコンピュータを用いて実現することもできる。図２０に示すような処理の手順を定めたコンピュータプログラムを記録した記録媒体を記録媒体読取装置で読み取り、読み取ったコンピュータプログラムをコンピュータに備えられたＲＡＭにロードし、コンピュータプログラムをＣＰＵ（プロセッサ）で実行することにより、コンピュータ上で脳状態推定装置５０を実現することができる。
　上述の実施の形態では、被験者の脳活動データとして、ＭＥＧデータを例に挙げて説明したが、脳活動データは、ＭＥＧデータに限定されるものではなく、ＥＥＧデータに対しても本実施の形態を適用することができる。
　脳状態推定装置５０は、生体磁気計測装置（ＭＥＧ計測装置）から、無線通信又は有線通信を経由して、ＭＥＧデータを取得することができる。ＭＥＧ計測装置は、ＭＥＧセンサ（検出コイル）を被験者の頭皮上に多数配置し、ＭＥＧセンサで記録された時系列データを出力する。多数のＭＥＧセンサそれぞれには、チャネルが対応付けられている。
　本実施の形態の脳状態推定装置は、時間軸を含む空間で変動する被験者の脳データを取得する取得部と、取得した脳データを変換して脳に関連する情報を生成する生成部と、生成した情報に基づいてエントロピーを算出する算出部と、算出したエントロピーがある確率測度を持つと仮定した場合のエントロピーを変換した際の測度の乖離度に基づいて前記被験者の脳状態を推定する推定部とを備える。
　本実施の形態の脳状態推定装置において、前記取得部は、脳活動データを含む脳データ、あるいは、脳活動データ、筋電位データ、心電図データ、眼電位データ、及び基準ノイズデータを含む脳データを取得し、前記生成部は、周波数変換又は非線形時系列解析による変換を用いて、律動の周波数帯域情報又はフラクタル情報を含む脳に関連する情報を生成し、前記算出部は、生成した律動の周波数帯域情報を含む情報に基づいて脳内の同期位相差量を特徴付けるエントロピー、又は生成したフラクタル情報を含む情報に基づいて時系列信号の乱雑性や位相構造の破れ具合を特徴付けるエントロピーを算出し、前記推定部は、算出したエントロピーの確率分布を平衡確率分布と仮定した場合のエントロピーの確率集合の乖離度に基づいて前記被験者の脳状態を推定する。
　本実施の形態の脳状態推定装置において、前記取得部は、脳内の複数の位置に対応する複数のチャネル毎の時系列の脳活動データを取得し、前記生成部は、前記複数のチャネルそれぞれの脳活動データ及び前記複数のチャネルのうちの２つのチャネル間の脳活動データに基づいて脳活動に関連する律動の周波数帯域情報を生成する。
　本実施の形態の脳状態推定装置において、前記算出部は、チャネル毎に、あるスペクトル密度が観測される測度に基づいてエントロピーを算出する。
　本実施の形態の脳状態推定装置において、前記算出部は、２つのチャネルで、あるスペクトル密度が同時に観測される測度に基づいてエントロピーを算出する。
　本実施の形態の脳状態推定装置において、前記推定部は、前記乖離度を健常者群に対する脳疾患群のマハラビノス距離で表現して脳状態を推定する。
　本実施の形態の脳状態推定装置は、前記取得部で取得した脳活動データから前記脳活動データの平均値を減算するバイアス補正、所定の周波数帯域のバンドパスフィルタ、所定の周波数のノッチフィルタ、所定期間毎の加算平均の少なくとも１つを用いて、前記脳活動データのノイズを除去する第１ノイズ除去部を備える。
　本実施の形態の脳状態推定装置は、前記取得部で取得した脳活動データに対して、独立成分分析及び信号空間投影法を用いて、前記脳活動データのノイズを除去する第２ノイズ除去部を備える。
　本実施の形態の脳状態推定装置は、脳内の複数の位置に対応する複数のチャネルのうちの２つのチャネル間のクロススペクトル分布を表示する表示処理部を備える。
　本実施の形態の脳状態推定装置において、前記表示処理部は、脳体積モデルに前記クロススペクトル分布を重畳して表示する。
　本実施の形態の脳状態推定装置において、前記表示処理部は、健常者と脳疾患を有する患者それぞれの前記クロススペクトル分布を表示する。
　本実施の形態の脳状態推定装置において、前記表示処理部は、健常者群に対する脳疾患群の脳状態の乖離度合いを表示する。
　本実施の形態のコンピュータプログラムは、コンピュータに、時間軸を含む空間で変動する被験者の脳データを取得し、取得した脳データを変換して脳に関連する情報を生成し、生成した情報に基づいてエントロピーを算出し、算出したエントロピーがある確率測度を持つと仮定した場合のエントロピーを変換した際の測度の乖離度に基づいて前記被験者の脳状態を推定する、処理を実行させる。
　本実施の形態の脳状態推定方法は、時間軸を含む空間で変動する被験者の脳データを取得し、取得した脳データを変換して脳に関連する情報を生成し、生成した情報に基づいてエントロピーを算出し、算出したエントロピーがある確率測度を持つと仮定した場合のエントロピーを変換した際の測度の乖離度に基づいて前記被験者の脳状態を推定する。
　本実施の形態のシステムは、時間軸を含む空間で変動する被験者の脳活動データを、脳機能の検査に係る多数の解を含む学習データを学習したモデルを用いて分類する機能を含む。
　本実施の形態のシステムは、時間軸を含む空間で変動する被験者の脳活動データから、脳の活動領域およびパターンを含む多数の入出力に係る解を学習したモデルを用いて、前記脳活動データにより示される信号源の時間軸を含む空間分布を推定する機能を含む。
　本実施の形態のシステムは、時間軸を含む空間で変動する被験者の脳活動データから推定された信号源の時間軸を含む空間分布と症例とを学習したモデルを用いて、被験者の脳状態を推定する機能を含む。
　本実施の形態のシステムにおいて、前記モデルは時系列ニューラルネットワークを含む。
　本実施の形態のシステムにおいて、前記脳活動データは、視覚、聴覚、触覚、味覚、嗅覚の少なくともいずれかを含む刺激に対する脳活動データを含む。
　本実施の形態の方法は、時間軸を含む空間で変動する被験者の脳活動データを、脳機能の検査に係る多数の解を含む学習データを学習した時系列モデルを用いて分類することを含む。
　本実施の形態の方法は、時間軸を含む空間で変動する被験者の脳活動データから、脳の活動領域およびパターンを含む多数の入出力に係る解を学習したモデルを用いて、前記脳活動データにより示される信号源の時間軸を含む空間分布を推定することを含む。
　本実施の形態の方法は、時間軸を含む空間で変動する被験者の脳活動データから推定された信号源の時間軸を含む空間分布と症例とを学習したモデルを用いて、被験者の脳状態を推定することを含む。
　本実施の形態のコンピュータプログラムは、前述のシステムとして稼働するための命令を有する。
　５０　脳状態推定装置
　５１　制御部
　５２　データ取得部
　５３　記憶部
　５４　前処理部
　５５　ＳＳＰ処理部
　５６　信号源位置推定部
　５７　信号分布算出部
　５８　スペクトル解析部
　５９　エントロピー算出部
　６０　脳状態推定部
　６１　表示パネル
　６２　操作部
　６３　表示処理部
 

Claims (23)

  1. 　時間軸を含む空間で変動する被験者の脳データを取得する取得部と、
    　取得した脳データを変換して脳に関連する情報を生成する生成部と、
    　生成した情報に基づいてエントロピーを算出する算出部と、
    　算出したエントロピーがある確率測度を持つと仮定した場合のエントロピーを変換した際の測度の乖離度に基づいて前記被験者の脳状態を推定する推定部と
    　を備える脳状態推定装置。
  2. 　前記取得部は、
    　脳活動データを含む脳データ、あるいは、脳活動データ、筋電位データ、心電図データ、眼電位データ、及び基準ノイズデータを含む脳データを取得し、
    　前記生成部は、
    　周波数変換又は非線形時系列解析による変換を用いて、律動の周波数帯域情報又はフラクタル情報を含む脳に関連する情報を生成し、
    　前記算出部は
    　生成した律動の周波数帯域情報を含む情報に基づいて脳内の同期位相差量を特徴付けるエントロピー、又は生成したフラクタル情報を含む情報に基づいて時系列信号の乱雑性や位相構造の破れ具合を特徴付けるエントロピーを算出し、
    　前記推定部は、
    　算出したエントロピーの確率分布を平衡確率分布と仮定した場合のエントロピーの確率集合の乖離度に基づいて前記被験者の脳状態を推定する、
    　請求項１に記載の脳状態推定装置。
  3. 　前記取得部は、
    　脳内の複数の位置に対応する複数のチャネル毎の時系列の脳活動データを取得し、
    　前記生成部は、
    　前記複数のチャネルそれぞれの脳活動データ及び前記複数のチャネルのうちの２つのチャネル間の脳活動データに基づいて脳活動に関連する律動の周波数帯域情報を生成する、
    　請求項２に記載の脳状態推定装置。
  4. 　前記算出部は、
    　チャネル毎に、あるスペクトル密度が観測される測度に基づいてエントロピーを算出する、
    　請求項３に記載の脳状態推定装置。
  5. 　前記算出部は、
    　２つのチャネルで、あるスペクトル密度が同時に観測される測度に基づいてエントロピーを算出する、
    　請求項３又は請求項４に記載の脳状態推定装置。
  6. 　前記推定部は、
    　前記乖離度を健常者群に対する脳疾患群のマハラビノス距離で表現して脳状態を推定する、
    　請求項２から請求項５のいずれか一項に記載の脳状態推定装置。
  7. 　前記取得部で取得した脳活動データから前記脳活動データの平均値を減算するバイアス補正、所定の周波数帯域のバンドパスフィルタ、所定の周波数のノッチフィルタ、所定期間毎の加算平均の少なくとも１つを用いて、前記脳活動データのノイズを除去する第１ノイズ除去部を備える、
    　請求項２から請求項６のいずれか一項に記載の脳状態推定装置。
  8. 　　前記取得部で取得した脳活動データに対して、独立成分分析及び信号空間投影法を用いて、前記脳活動データのノイズを除去する第２ノイズ除去部を備える、
    　請求項２から請求項７のいずれか一項に記載の脳状態推定装置。
  9. 　脳内の複数の位置に対応する複数のチャネルのうちの２つのチャネル間のクロススペクトル分布を表示する表示処理部を備える、
    　請求項１から請求項８のいずれか一項に記載の脳状態推定装置。
  10. 　前記表示処理部は、
    　脳体積モデルに前記クロススペクトル分布を重畳して表示する、
    　請求項９に記載の脳状態推定装置。
  11. 　前記表示処理部は、
    　健常者と脳疾患を有する患者それぞれの前記クロススペクトル分布を表示する、
    　請求項９又は請求項１０に記載の脳状態推定装置。
  12. 　前記表示処理部は、
    　健常者群に対する脳疾患群の脳状態の乖離度合いを表示する、
    　請求項９から請求項１１のいずれか一項に記載の脳状態推定装置。
  13. 　コンピュータに、
    　時間軸を含む空間で変動する被験者の脳データを取得し、
    　取得した脳データを変換して脳に関連する情報を生成し、
    　生成した情報に基づいてエントロピーを算出し、
    　算出したエントロピーがある確率測度を持つと仮定した場合のエントロピーを変換した際の測度の乖離度に基づいて前記被験者の脳状態を推定する、
    　処理を実行させるコンピュータプログラム。
  14. 　時間軸を含む空間で変動する被験者の脳データを取得し、
    　取得した脳データを変換して脳に関連する情報を生成し、
    　生成した情報に基づいてエントロピーを算出し、
    　算出したエントロピーがある確率測度を持つと仮定した場合のエントロピーを変換した際の測度の乖離度に基づいて前記被験者の脳状態を推定する、
    　脳状態推定方法。
  15. 　時間軸を含む空間で変動する被験者の脳活動データを、脳機能の検査に係る多数の解を含む学習データを学習したモデルを用いて分類する機能を含む、システム。
  16. 　時間軸を含む空間で変動する被験者の脳活動データから、脳の活動領域およびパターンを含む多数の入出力に係る解を学習したモデルを用いて、前記脳活動データにより示される信号源の時間軸を含む空間分布を推定する機能を含む、システム。
  17. 　時間軸を含む空間で変動する被験者の脳活動データから推定された信号源の時間軸を含む空間分布と症例とを学習したモデルを用いて、被験者の脳状態を推定する機能を含む、システム。
  18. 　請求項１５ないし１７のいずれかにおいて、
    　前記モデルは時系列ニューラルネットワークを含む、システム。
  19. 　請求項１５ないし１８のいずれかにおいて、
    　前記脳活動データは、視覚、聴覚、触覚、味覚、嗅覚の少なくともいずれかを含む刺激に対する脳活動データを含む、システム。
  20. 　時間軸を含む空間で変動する被験者の脳活動データを、脳機能の検査に係る多数の解を含む学習データを学習した時系列モデルを用いて分類することを含む、方法。
  21. 　時間軸を含む空間で変動する被験者の脳活動データから、脳の活動領域およびパターンを含む多数の入出力に係る解を学習したモデルを用いて、前記脳活動データにより示される信号源の時間軸を含む空間分布を推定することを含む、方法。
  22. 　時間軸を含む空間で変動する被験者の脳活動データから推定された信号源の時間軸を含む空間分布と症例とを学習したモデルを用いて、被験者の脳状態を推定することを含む、方法。
  23. 　コンピュータを請求項１５ないし１９のいずれかのシステムとして稼働するための命令を有するコンピュータプログラム。
PCT/JP2020/039091 2019-10-18 2020-10-16 脳状態推定装置、コンピュータプログラム、脳状態推定方法、脳機能の検査システムおよび方法 WO2021075548A1 (ja)

Priority Applications (1)

Application Number Priority Date Filing Date Title
JP2021552465A JPWO2021075548A1 (ja) 2019-10-18 2020-10-16

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2019191003 2019-10-18
JP2019-191003 2019-10-18

Publications (1)

Publication Number Publication Date
WO2021075548A1 true WO2021075548A1 (ja) 2021-04-22

Family

ID=75537413

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2020/039091 WO2021075548A1 (ja) 2019-10-18 2020-10-16 脳状態推定装置、コンピュータプログラム、脳状態推定方法、脳機能の検査システムおよび方法

Country Status (2)

Country Link
JP (1) JPWO2021075548A1 (ja)
WO (1) WO2021075548A1 (ja)

Cited By (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113288102A (zh) * 2021-06-11 2021-08-24 中国人民解放军陆军军医大学 一种无创监测脑血流的系统
CN113729735A (zh) * 2021-09-30 2021-12-03 上海交通大学 基于多域自适应图卷积神经网络的情感脑电特征表示方法
US11771377B1 (en) * 2022-12-28 2023-10-03 Neumarker, Inc Method and system for identifying cohorts of psychiatric disorder patients based on electroencephalograph

Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2008229238A (ja) * 2007-03-23 2008-10-02 National Institute Of Advanced Industrial & Technology 脳活動解析方法および装置
JP2009542351A (ja) * 2006-07-06 2009-12-03 リージェンツ オブ ザ ユニバーシティ オブ ミネソタ 時間的尺度を使用した脳パターンの解析
JP2015533559A (ja) * 2012-09-28 2015-11-26 ザ レジェンツ オブ ザ ユニヴァーシティー オブ カリフォルニア 知覚および認知プロファイリングのためのシステムおよび方法
US20180268942A1 (en) * 2016-12-06 2018-09-20 Darmiyan, Inc. Methods and systems for identifying brain disorders

Patent Citations (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2009542351A (ja) * 2006-07-06 2009-12-03 リージェンツ オブ ザ ユニバーシティ オブ ミネソタ 時間的尺度を使用した脳パターンの解析
JP2008229238A (ja) * 2007-03-23 2008-10-02 National Institute Of Advanced Industrial & Technology 脳活動解析方法および装置
JP2015533559A (ja) * 2012-09-28 2015-11-26 ザ レジェンツ オブ ザ ユニヴァーシティー オブ カリフォルニア 知覚および認知プロファイリングのためのシステムおよび方法
US20180268942A1 (en) * 2016-12-06 2018-09-20 Darmiyan, Inc. Methods and systems for identifying brain disorders

Cited By (4)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
CN113288102A (zh) * 2021-06-11 2021-08-24 中国人民解放军陆军军医大学 一种无创监测脑血流的系统
CN113729735A (zh) * 2021-09-30 2021-12-03 上海交通大学 基于多域自适应图卷积神经网络的情感脑电特征表示方法
CN113729735B (zh) * 2021-09-30 2022-05-17 上海交通大学 基于多域自适应图卷积神经网络的情感脑电特征表示方法
US11771377B1 (en) * 2022-12-28 2023-10-03 Neumarker, Inc Method and system for identifying cohorts of psychiatric disorder patients based on electroencephalograph

Also Published As

Publication number Publication date
JPWO2021075548A1 (ja) 2021-04-22

Similar Documents

Publication Publication Date Title
AU2019201766B2 (en) Systems and methods for diagnosis of depression and other medical conditions
Astolfi et al. Estimation of the effective and functional human cortical connectivity with structural equation modeling and directed transfer function applied to high-resolution EEG
WO2021075548A1 (ja) 脳状態推定装置、コンピュータプログラム、脳状態推定方法、脳機能の検査システムおよび方法
Astolfi et al. Assessing cortical functional connectivity by linear inverse estimation and directed transfer function: simulations and application to real data
Steele et al. Neuroimaging measures of error-processing: Extracting reliable signals from event-related potentials and functional magnetic resonance imaging
Samuelsson et al. Spatial fidelity of MEG/EEG source estimates: A general evaluation approach
Niso et al. What graph theory actually tells us about resting state interictal MEG epileptic activity
Kierkels et al. A model-based objective evaluation of eye movement correction in EEG recordings
Brookes et al. Task induced modulation of neural oscillations in electrophysiological brain networks
JP7221693B2 (ja) ヒトの脳の皮質機能の目録を作るための方法および磁気画像化デバイス
US20170332933A1 (en) Deep brain source imaging with m/eeg and anatomical mri
Schaworonkow et al. Is sensor space analysis good enough? Spatial patterns as a tool for assessing spatial mixing of EEG/MEG rhythms
Bardouille et al. Evidence for age-related changes in sensorimotor neuromagnetic responses during cued button pressing in a large open-access dataset
JPS63226340A (ja) 脳神経活動の位置と内部域の時間的関係を表示する方法とその装置
US8600493B2 (en) Method, apparatus and computer program product for automatic seizure monitoring
Mill et al. Empirical validation of directed functional connectivity
Zilber et al. Modulation of scale-free properties of brain activity in MEG
Sitnikova et al. Short timescale abnormalities in the states of spontaneous synchrony in the functional neural networks in Alzheimer's disease
US20160051161A1 (en) Method for locating a brain activity associated with a task
Hamandi et al. Non-invasive brain mapping in epilepsy: applications from magnetoencephalography
Barbati et al. Functional source separation from magnetoencephalographic signals
Khalid et al. Epileptic MEG spikes detection using amplitude thresholding and dynamic time warping
Vakorin et al. Complexity analysis of source activity underlying the neuromagnetic somatosensory steady-state response
Chiang et al. A generalized multivariate autoregressive (GmAR)-based approach for EEG source connectivity analysis
US11559244B2 (en) System and methods for consciousness evaluation in non-communicating subjects

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: 20875875

Country of ref document: EP

Kind code of ref document: A1

ENP Entry into the national phase

Ref document number: 2021552465

Country of ref document: JP

Kind code of ref document: A

NENP Non-entry into the national phase

Ref country code: DE

122 Ep: pct application non-entry in european phase

Ref document number: 20875875

Country of ref document: EP

Kind code of ref document: A1