System and method for pain monitoring using a multidimensional analysis of physiological signals

Abstract

The present invention is for a method and system for pain classification and monitoring optionally in a subject that is an awake, semi-awake or sedated.

Claims

1. A method for monitoring pain of a patient, the method comprising: obtaining at least two physiological signals comprising electrocardiogram (ECG) and at least one signal selected from the group consisting of: Photoplethysmograph (PPG), Galvanic Skin Response (GSR), Blood pressure, respiration, internal body temperature, skin temperature, electrooculography (EOG), pupil diameter, electroencephalogram (EEG), frontalis electromyogram (FEMG), electromyography (EMG), electro-gastro-gram (EGG), partial pressure of carbon dioxide, and accelerometer readings; processing said at least two physiological signals to improve signal quality, thereby obtaining at least two processed signals; generating a first vector, said first vector comprising at least three features extracted from said at least two physiological signals, wherein said at least three features comprise at least one feature selected from the group consisting of: RR/PQ/PR/QT/RS/ST interval, RR/PQ/PR/QT/RS/ST interval variability, Q/R/S/T/P amplitude, and any combination thereof; transforming said first vector into a second vector, said transformation comprising normalization; and monitoring the pain status of the patient by applying a classification algorithm adapted to classify said second vector into a graduated scale representing the level of pain; wherein said classification algorithm comprises an ensemble of classification and regression trees.

Description

BRIEF DESCRIPTION OF THE DRAWINGS

(1) The invention is herein described, by way of example only, with reference to the accompanying drawings. With specific reference now to the drawings in detail, it is stressed that the particulars shown are by way of example and for purposes of illustrative discussion of the preferred embodiments of the present invention only, and are presented in order to provide what is believed to be the most useful and readily understood description of the principles and conceptual aspects of the invention. In this regard, no attempt is made to show structural details of the invention in more detail than is necessary for a fundamental understanding of the invention, the description taken with the drawings making apparent to those skilled in the art how the several forms of the invention may be embodied in practice.

(2) In the drawings:

(3) FIG. 1 is schematic block diagram of an exemplary system according to the present invention.

(4) FIGS. 2A-C are schematic block diagram of the machine learning module of FIG. 1 in greater detail.

(5) FIGS. 3A-E are schematic block diagrams of optional embodiments of a system according to the present invention.

(6) FIG. 4 is an exemplary method according to the present invention.

(7) FIGS. 5A-C are scatter plots depicting pain classification according to an optional system and method of the present invention.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

(8) The present invention is of a system and a method for pain detection and monitoring most preferably detection and monitoring of pain is facilitated by processing physiological signals and features most preferably represented by a data vector comprising a great plurality of features. The principles and operation of the present invention may be better understood with reference to the drawings and the accompanying description.

(9) For the sake of clearly throughout the figures similar labels and numbering scheme is used throughout for equivalent or similarly functioning elements.

(10) FIG. 1 is a schematic block diagram of an exemplary system 100 according to the present invention for pain monitoring comprising physiological signal acquisition module 102 and processing module 104.

(11) Most preferably signal acquisition module 102 is provided for acquiring and measuring physiological measurements from a subject able to experience pain optionally including but not limited to a person is optionally provided with at least one or more preferably a plurality of sensor or transducers as are known and accepted in the art. Processing module 104 is most preferably provided with a plurality of sub-modules for example including but not limited to signal processing module 2, feature extraction module 3, machine learning module 4. Optionally and preferably signal processing module 2 provides for preprocessing of the acquired signal for example including but not limited to normalization, filtering, noise reduction, SNR optimization, domain transformations, statistical analysis, spectral analysis, wavelet analysis, or the like.

(12) Optionally and preferably feature extraction module 3 is provided for extracting features associated with the signals acquired with module 102. Optionally and preferably, feature extraction module may be further provided with an a priori data sub-module 6 most preferably for providing data beyond the acquired signals for example including but not limited to physician and/or caregiver data, medical history, genetic predispositions data, medical records or the like a priori data. Most preferably feature extraction is provided by performing signal processing techniques as is known and accepted in the art to extract from a physiological signal relevant and pertinent data. For example, a ECG measured acquired in module 102 is processed to provide a plurality of features. For example including but not limited to HRV, complex analysis, cardiac output data or the like. Most preferably feature extraction module 3 from the GPF vector for further analysis, for example including the features depicted in Table 1 above

(13) Machine learning module 4 is provided to classify the feature vector set provided from module 3. Optionally, machine learning module 4 may provide for feature selection and dimensionality reduction wherein the GPF vector undergoes feature selection and dimensionality reduction to provide a second vector therein providing a representation of the GPF in a smaller dimension. Optionally and preferably feature selection and dimensionality reduction techniques for example include but are not limited to Multi Dimensional Scaling (MDS), Principal Component Analysis (PCA), Sparse PCA (SPCA), Fisher Linear Discriminant Analysis (FLDA), Sparse FLDA (SFLDA), Kernel PCA (KPCA), ISOMAP, Locally Linear Embedding (LLE), Laplacian Eigenmaps, Diffusion Maps, Hessian Eigenmaps, Independent Component Analysis (ICA), Factor analysis (FA), Hierarchical Dimensionality Reduction (HDR), Sure Independence Screening (SIS), Fisher score ranks, t-test rank, Mann-Whitney U-test taken alone or in any combination thereof, or as known and accepted in the art.

(14) Optionally classification is performed on at least one of the GPF vector or a second vector comprising a dimensionally reduced version of GPF vector. Most preferably machine learning module 4 provides for pain classification into at least two or more classes, for example including pain and non-pain groups. Optionally, and preferably a plurality of optional classifiers may be used for example including but not limited Nearest Shrunken Centroids (NSC), Classification and Regression Trees (CART), ID3, C4.5, Multivariate Additive regression splines (MARS), Multiple additive regression trees (MART), Nearest Centroid (NC), Shrunken Centroid Regularized Linear Discriminate and Analysis (SCRLDA), Random Forest, Boosting, Bagging Classifier, AdaBoost, RealAdaBoost, LPBoost, TotalBoost, BrownBoost, MadaBoost, LogitBoost, GentleBoost, RobustBoost, Support Vector Machine (SVM), kernelized SVM, Linear classifier, Quadratic Discriminant Analysis (QDA) classifier, Naïve Bayes Classifier and Generalized Likelihood Ratio Test (GLRT) classifier with plug-in parametric or non-parametric class conditional density estimation, k-nearest neighbor, Radial Base Function (RBF) classifier, Multilayer Perceptron classifier, Bayesian Network (BN) classifier, multi-class classifier adapted from binary classifier with one-vs-one majority voting, one-vs-rest, Error Correcting Output Codes, hierarchical multi-class classification, Committee of classifiers or the like as are known and accepted in the art or any combination thereof.

(15) Optionally classification may be further provided with a graded scoring relating the level of pain classified.

(16) Output module 5 is provided to display or otherwise communicate the classification results provided by machine learning module 4. Optionally classification results may be displayed in a plurality of formats for example including printout, visual display cues, acoustic cues or the like. Optionally results may be displayed in graded or class formats, or the like.

(17) Optionally output module may communicate the classification results to and external device for further processing, medical intervention or the like, for example classification results may be communicated to a drug administration device to automatically, semi-automatically, or manually control the drug delivery of a pain medicament, or to indicate to a caregiver to control, change, decrease or increase the dosage or delivery of a pain reducing drug.

(18) FIGS. 2A-C provide a further optional depiction of machine learning module 4 of the processing module 104 of FIG. 1 wherein the feature vector set provided from module 3 is processed, most preferably sequentially, with a plurality of machine learning sub-modules, for example, preprocessing and normalization sub-module, feature selection and dimensionality reduction sub-module and classification sub-module. Machine-learning sub-modules implements machine learning techniques and methods as is known and accepted in the art. Optionally, sub-modules of the machine learning module comprising preprocessing and normalization, feature selection and dimensionality reduction may be reorganized and provided in all possible combination.

(19) FIG. 2B provides a depiction of a non limiting example of an optional embodiment of FIG. 2A, also described in Example 1 below, wherein the machine learning module 4 of the processing module 104 of FIG. 1 is provided by a collection of machine learning sub-modules. For example, preprocessing is provided by a plurality of actions including outlier removal and signal normalization utilizing baseline zero mean unit variance; feature selection is provided through the use of a Fisher score; dimensionality reduction is provided by Sparse Fisher Linear Discriminant Analysis (SFLDA) while classification is provided by a RealAdaboost classifier. Most preferably machine learning module 4 of FIG. 1 comprises at least one preprocessing sub-module, at least one feature selection and dimensionality reduction sub-module, and at least one classification sub-module; of the optional techniques described above, in any combination thereof as is know and accepted in the art.

(20) FIG. 2C provides a depiction of a non limiting example of an optional embodiment of FIG. 2A, also described in Example 3 below, wherein the machine learning module 4 of the processing module 104 of FIG. 1 is provided by a collection of machine learning sub-modules. For example, preprocessing is provided by a plurality of actions including outlier removal and signal normalization utilizing baseline zero mean unit variance; feature selection and dimensionality reduction is provided by Hierarchical dimensionality reduction while classification is provided by a Random Forest classifier.

(21) FIGS. 3A-E depicts optional hardware configurations of optional systems according to optional embodiments of the present invention as described in FIG. 1 above. FIGS. 3A-3E depicts different remote and local configurations of signal acquisition, processing and display. FIGS. 3A and 3E depicts systems that may be adapted to provide for a fully local system, FIGS. 3B and 3C depicts systems that may be optionally configured to be semi remote, while FIG. 3D depicts a system that may be adapted to be fully remote.

(22) FIG. 3A depicts optional configuration system 300 comprising computer 302 and an external signal acquisition and display monitor 308, optionally and preferably communicating using communication protocol 310. Most preferably, monitor 308 comprises a stand alone or external monitor comprising signal acquisition and display that optionally and preferably provides for signal acquisition that is preferably processed with computer 302. Optionally and preferably computer 302 provides for the processing for pain monitoring and detection according to the present. invention. Optionally computer 302 may be realized as a processor, server or the like computing device as is known and accepted in the art. Optionally communication protocol 310 mediates communication and data exchange between computer 302 and external monitor 308 providing for pain classification. Optionally, communication protocol may be provided in a plurality of optional communication protocols as is known and accepted in the art for example including but not limited to wireless, wired, cellular, internet, Bluetooth, optical, IR or the like communication protocols as is known and accepted in the art. Optionally and preferably computer 302 provides for pain classification and monitoring according to the present invention while physiological signal acquisition and display is provided by monitor 308 most preferably monitor 308 acquires the signals, transmit them to processor 302 and optionally and preferably displays the classification result

(23) FIG. 3C depicts an optional system 303 similar in configuration to system 300 of FIG. 3A above wherein the physiological signals are preferably provided by both monitor 308 and acquisition transducers 304 wherein both are communicated to computer 302 for processing according to an optional method according to the present invention.

(24) FIG. 3B depicts optional configuration system 301 comprising computer 302, signal acquisition transducers 304, physiological signal acquisition device 306, acquiring and displaying monitor 308 and communication protocol 310. Optionally signal acquisition transducer 304 is most preferably provided in the form of optional sensors dependent on the signal being acquired that is preferably coupled to signal acquisition device 306 most preferably to sample, amplify and process and record a physiological signal. Optionally and preferably device 306 is a mobile, optionally provided in a plurality of formats for example including but not limited to mobile telephone, PDA, hand held device, MP3 player, dedicated or converted device adept for receiving and communicating a physiological signal. Optionally device 306 comprises a processor for performing initial signal processing of the physiological signal. Optionally, device 306 does not comprise a processor adept for processing the sampled physiological signals and is linked using communication protocol 310 to higher processing centers, for example computer 302. Optionally, device 306 and computer 302 comprise a master slave processing protocols for processing the physiological signals sampled with system acquisition transducers 304. Optionally communication protocol 310 may be further facilitated through an interne connection. Optionally, computer 302 and display 308 may remote to the signal acquisition transducer 304 and device 306 relaying on communication protocol 310 to communicate therebetween. Optionally, computer 302 and acquiring and/or displaying monitor 308 may be implemented as a call center, telemedicine center, emergency medical center, or remote center for remote pain classification and detection for pain management.

(25) FIG. 3D depicts system 305 for remote pain classification comprising a remote signal acquisition device 306 is utilized to acquire and communicate physiological signals for pain classification provided for with computer 302 provided in the form of a PDA comprising an intrinsic display. System 305 therefore provides a remote system comprising remote signal acquisition as well as remote signal processing and display. Conversely, FIG. 3E depicts system 307 comprising both local signal acquisition and processing.

(26) FIG. 4 shows a flowchart of an exemplary method according to the present invention for pain monitoring and classification. In stage 1 physiological signals are acquired, preferably a plurality of signals are obtained as previously described in FIG. 1 using a plurality of optional sensors and/or transducers. Next in stage 2, signal preprocessing is performed for example including but not limited to SNR optimization, signal normalization, filtering or the like. Next in stage 3, features extraction is performed to abstract and to preferably derive a great plurality of features GPF from the acquired physiological signals. Next classification processed is initiated with a training process on the training set as depicted in stages 4′ and 5′. Most preferably an initial training process is performed in stage 4′ wherein feature selection and dimensionality reduction functions are trained to implement the optional dimensionality reduction and feature selection techniques as previously described. Next in stage 5′ the classifier is trained and implemented to identify the pain classification of the training set.

(27) Optionally and preferably once the classifier has been trained and set following stage 5′ future features obtained from the physiological signals are classified without further training of the classifying system of the present invention and implemented in stages 4 and 5. In stage 4 feature selection and dimensionality reduction is performed according to the dimensionality training determined in stage 4′. Next in stage 5 classifications is performed with a classifier to determine the pain according to the dimensionality reduced GPF vector. Finally in stage 6 pain monitoring and classification is displayed according to the appropriate time scale.

EXAMPLE 1: PAIN CLASSIFICATION UTILIZING A COMBINATION OF FISHER SCORE+SFLDA+REALADABOOST

(28) The following is a non limiting example of the implementation of pain monitoring according to an optional embodiment of an optional method for pain monitoring and classification according to the present invention. The following example provides an illustrative example of pain classification for pain monitoring in a subject wherein the processing methods comprise utilizing a priori data, feature extraction and selection based on Fisher Score rank, dimensionality reduction using SFLDA, and classification using a RealAdaboost classifier within a Boosting framework.

Materials and Methods

(29) An experiment was conducted in order to develop, validate performances and evaluate efficacy of the classification method and system for non-invasive automated pain monitoring according to the present invention. Primary outcome measures was to compare the pain monitoring results with the subjective pain report measured by the visual Numeric Pain Scale (NPS) to a given pain stimulus with reports divided to pain/no-pain (binary test).

(30) The study included 26 healthy volunteers. Inclusion criteria were: (I) free from chronic pain of any type, (II) no medication usage except for oral contraceptives, (III) ability to understand the purpose and instructions of the study, and (IV) blood pressure<140/90. Exclusion criteria were: (I) any type of preexisting condition, (II) use of medications or recreational drugs, or (III) pregnancy. The study was approved by the local ethics committee, and a written informed consent was obtained from all participants prior to the beginning of the experiment.

(31) The cold pressor test (CPT) and heat pain test were chosen to provoke the experimental pain. The cold pressor test apparatus (ChillSafe 8-30, ScanLaf A/S Denmark) is a temperature-controlled water bath with a maximum temperature variance of ±0.5° C., which is continuously stirred by a pump. Volunteers were asked to place their right foot (until above the ankle) in the CPT bath in a still position and maintain their foot in the water for 1 min each session. A thermal testing analyzer (TSA) thermode of 30×30 mm (Medoc TSA-2001 device, Ramat Ishai, Israel) was attached to the skin of the right forearm to initiate heat pain. During pain stimuli sessions the thermode was heated at 10° C./sec to the target temperature (39° -48.5° C.) and with a plateau lasting 60 sec.

(32) In order to evaluate the volunteer's cold and heat pain sensitivity they were exposed to a range of different temperatures and reported perceived pain on a 0-100 numeric pain scale (NPS). For evaluation of subject's heat pain sensitivity they were exposed to 11 heat stimuli, ranging from 37° C. to 50° C. with increasing rate of 10° C./sec, each with a plateau lasting for 10 seconds. In order to evaluate subject's cold pain sensitivity we expose subjects for one min to the CPT using water temperature of 12° C. Subjects reported their pain every 10 sec. The cold and heat apparatus were then appropriately calibrated to initiate feeling of no-pain (NPS 0-30), and pain (NPS>70). Due to the limitations of the GCP and ethics committee approval, the minimum/maximum temperatures used were 1° C. for the cold pain and 48.5° C. for the heat pain.

(33) Subjects received a full explanation about the purpose and design of the study and signed a written informed consent form. Prior to the beginning of the experiment, the familiarizing and calibration sessions were conducted. Next, the experimenter connected the sensors to the subject and during the next 5 min physiological signals were rerecorded for baseline normalization purposes. During the two sessions each participant received 4 heat stimuli and 3 cold stimuli sessions, lasting for 1 minute each with intervals of 10-15 minutes between stimuli and 30-45 minutes between sessions. Volunteers were unaware of the stimuli intensity. The order of the heat pain stimuli was randomly assigned, while the order of the cold stimuli was progressive from low to high level of pain (to avoid adaptation). The “no-pain” stimuli (25° C. CPT and 39° C. heat sensor) was introduced with intention to stimulate a sensory experience similar to pain sessions but without painful stimulus.

(34) The NPS report, the physiological signals and the intensity of stimuli were recorded and synchronized for subsequent processing.

(35) The physiological signals were recorded and stored in a personal computer by the BioPac MP 100 system (BioPac System Inc., Calif., USA) and its companion software AcqKnowledge 3.9.1 (BioPac System Inc.). One-lead electrocardiogram (ECG) signal, 2-channel electroencephalogram (EEG) signal from forehead, photoplethysmograph (PPG) signal from right hand finger, and one lead external electromyogram (EMG) signal from right trapezius muscle were sampled with a frequency of 500 Hz. External skin temperature from the dorsum of the right hand, respiration, and galvanic skin response (GSR) from right hand fingers were sampled with a frequency of 32.5 Hz. In addition, continuous blood pressure signal was non-invasively measured using the Finometer MIDI (Finapres Medical Systems BV, Amsterdam, The Netherlands) and data were recorded using companion software BeatScope EASY (Finapres Medical Systems BV). Continuous blood pressure was sampled with a frequency of 200 Hz.

(36) The recorded physiological signals were extracted, synchronized and processed in off-line way using Matlab® 2009 scientific software (The Mathworks, Inc., Mass., USA).

(37) The data was processed offline on PC computer. All signals were processed using routine signal processing methods for noise and artifact filtering (Oppenhem & Shafer 1999). For some signals (EEG, ECG, etc.) were used signal-specific data processing methods (Rangayyan 2002, Sannei & Chambers 2007). All extracted parameters were averaged (if applicable) with non-overlapping windows of 10 sec. Features utilized are summarized in Table 1 above were extracted from a plurality of raw physiological data. The extracted features were deployed to training and validate classification algorithm.

(38) The machine learning module implemented as described in greater detail in FIG. 2B above.

(39) During the preprocessing and normalization stages the data was manually examined by trained professional. Patients with over-sensitivity and under-sensitivity were removed from training data. This exclusion from training data set is a non-limiting example of a priori knowledge that is incorporated into the system and method of the present invention most preferably to improve classification results. Nevertheless, the data excluded from training the optional classifier was used to test the classifier.

(40) Next, feature normalization was performed by removing the patient's features baseline mean and normalizing the patient's feature baseline variability, in the following manner:

(41) $\frac{X_i-\mathrm{av}g\left(X_i^{\mathrm{baseline}}\right)}{std\left(X_i^{\mathrm{baseline}}\right)}$

(42) Next Feature Selection step was performed. During the training phase features with maximum Fisher score were selected. Fisher score of i'th feature was defined as

(43) $F_i=\frac{{\left(\mathrm{avg}\left(X_j^{\mathrm{pain}}\right)-\mathrm{avg}\left(X_i^{\mathrm{no}\text{-}\mathrm{pain}}\right)\right)}^2}{\mathrm{var}\left(X_i^{\mathrm{pain}}\right)+\mathrm{var}\left(X_i^{\mathrm{no}\text{-}\mathrm{pain}}\right)}$

(44) The 100 best features with highest Fisher scores were predetermined during the training phase. During feature selection only the 100 highly ranked features were kept, the remaining features were removed from consideration.

(45) Next, Dimensionality Reduction was performed. During the training phase of the non-limiting and optional method according to the present invention a transformation matrix was determined based on the principle of Sparse Fisher Linear Discriminant Analysis (SFLDA). The methods for calculation of a single SFLDA transformation vector were rigorously described in (Moghaddam, et al., 2006). wherein the transformation vector is calculated either by exhaustive search or using greedy algorithms (forward or backward). Multiplication of the transformation vector by a vector of the selected features following feature selection as described above, results in a number which represents the first reduced discriminate dimension. For pain monitoring unique specification, novel iterative procedure for finding multiple transformation vectors and discriminative components was used as follow: Input: B—N×N pooled between-class covariance matrix, W—N×N pooled within-class covariance matrix, k—required sparsity level for transformation vectors, P—number of transformation vectors Algorithm: 1. Set p=1, l.sup.p={1 . . . N} 2. Input B.sub.l.sub.p, W.sub.l.sub.p into SFLDA algorithm in order to compute optimum sparsity pattern I.sub.opt.sup.p sparse transformation vector α.sub.opt.sup.p. 3. Set l.sup.p+1=l.sup.p/l.sub.opt.sup.p. Set p=p+1. If p<P go back to step 2. Output: Set of sparsity patterns {I.sub.opt.sup.1, . . . , I.sub.opt.sup.P} and transformation matrix A=[α.sub.opt.sup.1 . . . α.sub.opt.sup.P].

(46) Transformation matrix A was multiplied by a vector of selected features and transformed it into a new vector most preferably having a significantly reduced dimension. This modification of SFLDA is known as Sparse Linear Discriminant Component Analysis (SLDCA). In an experiment each transformation vector had sparsity level 5, and a total 10 such vectors were computed. Therein according to this non-limiting example a feature set of 100 was reduced to 10 following an optional dimensionality reducing step comprising SFLDA as described above

(47) The final stage of classification according to an optional method of the present invention was set to classify the vector of reduced dimensions into a binary class of Pain and No-Pain classes. An optional classifier according to the present invention, RealAdaboost was chosen and trained during a training phase.

(48) In order to assess performances of proposed method and apparatus the leave-one-out cross-validation (Hastie et al., 2009) scheme has been used. In order to prevent the situation where the algorithm was both trained and validated on the same data (in our case same subject), the algorithm was applied N times, where N is the number of subjects. At each run, data were included in the training set from all subjects excluding one, and then the trained algorithm was scored on data from this subject. The Test Error is estimated by averaging over classification errors from each of the N runs.

(49) The algorithm was also tested on patients, which were declared as outliers. In such cases, the algorithm was trained on all non-outlier patients and was tested on all outliers.

(50) Results

(51) The performance of the algorithm is presented in the Table 2. The overall agreement, sensitivity, and positive predictive values (PPV) are presented, together with their respective 95% exact binomial confidence intervals

(52) TABLE-US-00002 TABLE 2 Results Pain level estimated pain no pain Total Pain level pain 425 68 493 no pain 181 1528 1709 Total 606 1596

(53) TABLE-US-00003 TABLE 3 Pain vs. No pain Percent 95% exact binomial confidence interval Overall agreement 88.69% 87.29% 89.99% PPV pain 70.13% 66.31% 73.75% PPV no pain 95.74% 94.63% 96.68% Sensitivity pain 86.21% 82.84% 89.13% Sensitivity no pain 89.41% 87.85% 90.83%

EXAMPLE 2: FEATURE SELECTION AND DIMENSIONALITY REDUCTION

(54) The following example is of a non-limiting example showing the importance of machine learning techniques according to the present invention to improve and provide for pain monitoring, classification and identification. During the clinical study described in Example 1 above all features described in Table 1 were measured. Physiological signals such as an accelerometer, SPO2 were excluded as the clinical study was performed on awake patients. Fisher scores were calculated to rank all available features, according to an optional method of the present invention.

(55) Analysis of the features showing the following features as having a high Fisher score, as summarized in Table 4 and FIGS. 5A-C.

(56) As can be seen these parameters represent both an autonomic system activity of a patient (e.g. amplitudes of PPG signal, Count of Spontaneous Fluctuations in GSR signal, H.R.V extracted from PPG or ECG and its spectral analysis) as well as behavioral activity, such as respiration rate and deviation, muscle activity (which may suggest on a certain discomfort).

(57) TABLE-US-00004 TABLE 4 Total Signal Temporal Feature Spectral Features Features GSR Count of Spontaneous Power of frequency 8 Fluctuations (Peaks), bands VLF/ Weighed Sum of Peaks LF/HF/above HF (weighted by their amplitude), Signal Max value in window, Smoothed Signal Value in window PPG Mean and Standard Power of frequency 9 Deviation of amplitude bands HF and of Peak, Trough, and LF, and LF/HF Max Rate point. Mean ratio calculated and Standard Deviation from spectral of AUC, Peak to Notch analysis of Peak interval, and Pulse to Notch Transition Time (Peak interval and Pulse to R-ECG interval). Transition Time. Respiratory Mean and Standard Power of frequency Deviation of the bands VLF/ amplitude of Peak, Peak LF/HF/above HF to Peak interval (respiratory rate), Smoothed Signal value, Signal Max and Min values in Window Temperature Weighed Sum of Peaks Power of frequency 9 (fluctuations) in bands VLF/ window, Mean LF/above HF. amplitude of Peaks in window ECG Mean and Standard Power of HF and LF 4 Deviation of R peak bands from Heart rate variability EEG/EMG Barlow's mean Energy in band 9 amplitude, Hjorth's 30-70 Hz, mobility, between mean power, channels covariance (we Barlow's mean use 2-channel EEG), frequency and time differentials of a spectral purity signal. EGG Power of below 1 LF band

(58) For example, FIG. 5A depicted scatter plots of 5 random features showing two classes of pain namely pain (“o” open circles and) and non pain (“+” plus sign) for a classification attempted with the initial features that have not undergone dimensionality reduction according to the present invention. Accordingly, inspection of the scatter plots clearly shows that the two classes are not separable either by linearly or non-linearly classification.

(59) FIG. 5B depicts scatter plots of 5 features extracted from the GPF vector based on the highest Fisher scores and as depicted in Table 4: 1st Fisher feature is (combination): number of Peaks, Weighed Sum of Peaks, Signal Max value in window, Smoothed Signal Value in window of GSR signal. 2nd Fisher feature is: Power of HF band of GSR signal. 3rd Fisher feature is: Standard deviation of Peaks to peak intervals of Respiratory signal. 4th Fisher feature is combination: Mean of Peak, Maximum Rate Point, and Trough Heart Rate from PPG signal. 5th Fisher feature is: Standard deviation of amplitude of Peaks of Respiratory signal.

(60) Although inspection of FIG. 5B shows that scoring features provides improved classification as can be seen by the reduced overlap between the two classes depicted when compared to FIG. 5A however a substantial degree of overlap is observed.

(61) FIG. 5C depicts scatter plots of 5 SLDCA components wherein dimensionality reduction technique utilized SLDCA substantially improves pain classification between the two classes.

(62) 1st SLDCA component comprises a combination of: Standard deviation of Peaks to peak intervals of Respiratory signal (3 Fisher feature) Deviation of amplitude of Peaks of Respiratory signal (5 Fisher feature) Number of Peaks, Weighed Sum of Peaks, Signal Max value in window, Smoothed Signal Value in window of GSR signal (1 Fisher feature) Power of LF band of PPG signal Mean of Peak, Maximum Rate Point, and Trough, Heart Rate from PPG signal (4 Fisher feature).

(63) 2nd SLDCA component comprises a combination of: Mean of the amplitude of Peak of Respiratory signal Standard Deviation of amplitude of Peak of Respiratory signal Power of bellow LF band EGG Between channels covariance EEG Power of above HF band of GSR

(64) 3rd SLDCA component is a combination of: First differential of EMG signal Standard Deviation of Peak-to-Peak interval of EGG Peak-to-Notch interval of PPG. Power of HF band of peak-to-notch variability. Power of VLF band Temperature

(65) 4th SLDCA component is a combination of: Power of HF band of Respiratory signal. Heart Rate from ECG signal Standard Deviation of amplitude of Peak of EGG Pulse transition time (Peak to R) from ECG and PPG Standard Deviation of Peak-to-Notch interval of PPG

(66) 5th SLDCA component is a combination of: Power of bellow VLF band of Respiratory Power of LF band EGG Standard Deviation of Pulse transition time (Peak to R) from ECG and PPG LF/HF ratio of PPG Power of LF band Temperature

(67) Although inspection of FIG. 5C shows that SLDCA dimensionality reduction provides improved classification as can be seen by the reduced overlap between the two classes depicted when compared to FIG. 5B still a substantial degree of overlap is observed, therefore classification module is still required.

EXAMPLE 3: PAIN CLASSIFICATION UTILIZING A COMBINATION OF HIERARCHICAL DIMENSIONALITY REDUCTION+RANDOM FOREST CLASSIFICATION

(68) The following example is of a non-limiting example showing the importance of machine learning techniques according to the present invention to improve and provide for pain monitoring, classification and identification, as depicted in FIG. 2C wherein the system and method of the present invention comprises a wherein feature selection and dimensionality reduction is provided by Hierarchical Dimensionality Reduction (HDR), and classification is provided by a Random Forest classifier.

Subjects

(69) The study included 36 healthy volunteers, 23 female and 13 male, aged 20 to 38 years (mean (SD) 26(4.3)). Inclusion criteria were: (I) free from chronic pain of any type, (II) no medication usage except for oral contraceptives, (III) ability to understand the purpose and instructions of the study, and (IV) blood pressure<140/90. Exclusion criteria were: (I) any type of preexisting condition, (II) use of medications or recreational drugs, or (III) pregnancy. The study was approved by the local ethics committee, and a written informed consent was obtained from all participants prior to the beginning of the experiment.

Instruments for Pain and Stress Stimulation

Assessment of Cold and Heat Pain Perception

(70) The cold pressor test (CPT) and heat pain test were chosen to provoke the experimental pain. The cold pressor test apparatus (ChillSafe 8-30, ScanLaf A/S Denmark) is a temperature-controlled water bath with a maximum temperature variance of ±0.5° C., which is continuously stirred by a pump. Volunteers were asked to place their right foot (until above the ankle) in the CPT bath in a still position and maintain their foot in the water for 1 min each session. A thermal testing analyzer (TSA) thermode of 30×30 mm (Medoc TSA-2001 device, Ramat Ishai, Israel) was attached to the skin of the right forearm to initiate heat pain.

(71) During pain stimuli sessions the thermode was heated at 10° C./sec to the target temperature (39°-48.5° C.) and with a plateau lasting 60 sec.

(72) In order to evaluate the volunteer's cold and heat pain sensitivity they were exposed to a range of different temperatures and reported perceived pain on a 0-100 numeric pain scale (NPS). For evaluation of subject's heat pain sensitivity they were exposed to 11 heat stimuli, ranging from 37° C. to 50° C. with increasing rate of 10° C./sec, each with a plateau lasting for 10 seconds. In order to evaluate subject's cold pain sensitivity we expose subjects for one min to the CPT using water temperature of 12° C. Subjects reported their pain every 10 sec. The cold and heat apparatus were then appropriately calibrated to initiate feeling of no-pain (NPS 0-15), mild pain (NPS 15-45), moderate pain (NPS 45-75), and severe pain (NPS>75). Due to the limitations of the GCP and ethics committee approval, the minimum/maximum temperatures used were 1° C. for the cold pain and 48.5° C. for the heat pain. On several occurrences the calibrated temperatures of “moderate pain” were significantly close to minimum/maximum allowed temperatures. At such cases “severe pain” stimuli temperatures were set to minimum/maximum temperatures.

Mental Stress Protocol

(73) The control sessions of mental stress was performed. In order to induce mental stress the mental arithmetic test Paced Auditory Serial Addition Task (PASAT), (Gronwall & Wrightson 1974) has been used, which consists of 1 min of auditory presentation of random digits from one to nine, with an interval of 2 sec. The subjects' were asked continuously express the sum of the two last digits.

Physiological Signals Recording

(74) The physiological signals were recorded and stored in a personal computer by the BioPac MP 100 system (BioPac System Inc., Calif., USA) and its companion software AcqKnowledge 3.9.1 (BioPac System Inc.). One-lead electrocardiogram (ECG) signal, 2-channel electroencephalogram (EEG) signal from forehead, photoplethysmograph (PPG) signal from right hand finger, and one lead external electromyogram (EMG) signal from right trapezius muscle were sampled with a frequency of 500 Hz. External skin temperature from the dorsum of the right hand, respiration, and galvanic skin response (GSR) from right hand fingers were sampled with a frequency of 32.5 Hz. In addition, continuous blood pressure signal was non-invasively measured using the Finometer MIDI (Finapres Medical Systems BV, Amsterdam, The Netherlands) and data were recorded using companion software BeatScope EASY (Finapres Medical Systems BV). Continuous blood pressure was sampled with a frequency of 200 Hz.

Data Processing and Parameters Extraction

(75) The recorded physiological signals were extracted, synchronized and processed in off-line way using Matlab® 2009 scientific software (The Mathworks, Inc., Mass., USA).

(76) All signals were processed using routine signal processing methods for noise and artifact filtering (Oppenhem & Shafer 1999). For some signals (EEG, ECG, etc.) were used signal-specific data processing methods (Rangayyan 2002, Sannei & Chambers 2007). All extracted parameters were averaged (if applicable) with non-overlapping windows of 10 sec. The detailed list of parameters which were extracted from above mentioned physiological signals can be found in Table 1

Normalization and Averaging

(77) All continuous parameters were normalized by removing the parameter's baseline mean and normalizing the parameter's baseline variability, in the following manner:

(78) 0$\frac{X_i-\mathrm{av}g\left(X_i^{\mathrm{baseline}}\right)}{std\left(X_i^{\mathrm{baseline}}\right)}$

(79) In order to avoid bias and over-fitting, baseline signals of a specific subject was not used for normalization of the subject's parameters. In other words parameters of the subject where normalized based on baseline parameters of all other volunteers. Categorical parameters were not normalized.

(80) Finally, the parameters during each stimulus were averaged. Thus each 1 min stimulus was represented by a vector of normalized parameters extracted from recorded physiological signals.

Dimensionality Reduction

(81) All parameters with correlation r>0.8 were identified and consequently averaged by Hierarchical Dimensionality Reduction (HDR) method Duda et al. 2000.

Classification Algorithm: Random Forest

(82) The goal of the presented example was to demonstrate abilities of a pain monitoring device which optionally and preferably provides for predicting the presence of pain sensation and classify its level. A first vector, according to optional embodiments of the present invention, comprising temporal, spectral, and other parameters extracted during data processing and parameter extraction step constitutes prediction variables which we will use in prediction procedure.

(83) In this example the classification algorithm which has been chosen was the Random Forest. Random Forest (Breiman 2001b) is a statistical learning procedure that makes a prediction by aggregating the results from an ensemble of classification and regression trees (CART) (Breiman et al. 1984). Random Forest averages over multiple CART trees increase the stability of the final algorithm. Each tree is grown on a bootstrapped sample (random sampling with repetition) of original training data. During the growing process of each tree, features are randomly sampled as well to find a best split of each node of the tree.

Classification Algorithm Testing Methodology

(84) In order to assess performances of proposed method and apparatus the leave-one-out cross-validation (Hastie et al., 2009) scheme has been used. In order to prevent the situation where the algorithm was both trained and validated on the same data (in our case same subject), the algorithm was applied N times, where N is the number of subjects. At each run, data were included in the training set from all subjects excluding one, and then the trained algorithm was scored on data from this subject. The Test Error is estimated by averaging over classification errors from each of the N runs.

Study Design

(85) Subjects received a full explanation about the purpose and design of the study and signed a written informed consent form. Prior to the beginning of the experiment, the familiarizing and calibration sessions were conducted. Next, the experimenter connected the sensors to the subject and during the next 5 min physiological signals were rerecorded for baseline normalization purposes. After the baseline the physical stress session was performed. During the two sessions each participant received 3 heat stimuli and 3 cold stimuli sessions (mild pain, moderate pain and severe pain), lasting for 1 minute each with intervals of 10-15 minutes between stimuli and 30-45 minutes between sessions. Volunteers were unaware of the stimuli intensity. The order of the heat pain stimuli was randomly assigned, while the order of the cold stimuli was progressive from low to high level of pain (to avoid adaptation). The “no-pain” stimuli (25° C. CPT and 39° C. heat sensor) was introduced with intention to stimulate a sensory experience similar to pain sessions but without painful stimulus.

Results

(86) The performance of the algorithm is assessed by presenting tables of the stimuli inflicted on the volunteers by experimenter versus the pain as predicted and classified by the classification software based solely on recorded physiological signals.

(87) Since the cost function was used in the training of the classification algorithm, the weight matrices should be used to rank the estimation error, i.e., to give a smaller penalty to a small error than a large one; for example, in case of high pain stimulation, the device estimation of “no pain” will be considered a worse error than a “mild pain” estimation. If pain levels are ranked by their severity (increasing or decreasing order) the “weight” of estimated pain level i given stimulated pain level j are defined according to the formula:
w(i,j)=1−abs(i−j)/(N−1),
where N is a total number of pain levels, and abs( ) denotes absolute value. This formula intends to introduce ranking constraint into pain level estimation.

(88) The weighting matrix used in the case of three categories is:

(89) TABLE-US-00005 Pain level estimated Pain level stimulated Severe Mild No Pain Severe 1 1/2 0 Mild ½ 1 1/2 No Pain 0 1/2 1
In case of two categories no weights are used, in other words all errors had equal costs. In these cases the overall agreement, sensitivity, and positive predictive values (PPV) are presented, together with their respective 95% exact binomial confidence intervals. Weights are introduced in sensitivity and PPV calculations for multiclass cases. Thus, the resulted sensitivity and PPV values along with their corresponded CI values are weighted sensitivity, PPV and CI respectively. In these cases, the weighted overall agreement, weighted sensitivity and weighted PPV are presented, together with their respective 95% exact Wilson-score confidence intervals.

(90) TABLE-US-00006 TABLE 5 Severe pain vs. No pain Pain level estimated severe pain no pain Total Pain level inflicted severe pain 249 33 282 no pain 28 271 299 Total 277 304 581

(91) TABLE-US-00007 TABLE 6 Severe pain vs. No pain 95% exact binomial Percent confidence interval Overall agreement 89.50% 86.72% 91.87% PPV severe pain 89.89% 85.72% 93.18% PPV no pain 89.14% 85.09% 92.41% Sensitivity severe pain 88.30% 83.96% 91.81% Sensitivity no pain 90.64% 86.75% 93.69%

(92) TABLE-US-00008 TABLE 7 Severe pain vs. Mental Stress Pain level estimated severe pain mental stress Total Pain level inflicted severe pain 292 37 329 mental stress 69 260 329 Total 361 297 658

(93) TABLE-US-00009 TABLE 8 Severe pain vs. Mental Stress 95% exact binomial Percent confidence interval Overall agreement 83.89% 80.85% 86.62% PPV severe pain 80.89% 76.44% 84.81% PPV mental stress 87.54% 83.24% 91.07% Sensitivity severe pain 88.75% 84.83% 91.96% Sensitivity mental stress 79.03% 74.22% 83.30%

(94) TABLE-US-00010 TABLE 9 Severe pain vs. Mild Pain vs. No pain Pain level estimated severe pain mild pain no pain Total Pain level inflicted severe pain 191 75 16 282 mild pain 59 178 39 276 no pain 4 69 226 299 Total 254 322 281 857

(95) TABLE-US-00011 TABLE 10 Severe pain vs. Mild Pain vs. No pain 95% exact Wilson Percent confidence interval Overall agreement 83.55% 80.92% 85.88% PPV severe pain 86.81% 82.10% 90.43% PPV mild pain 77.64% 72.78% 81.85% PPV no pain 87.37% 82.97% 90.75% Sensitivity severe pain 81.03% 76.05% 85.18% Sensitivity mild pain 82.25% 77.30% 86.30% Sensitivity no pain 87.12% 82.85% 90.45%

(96) These results show that the system and method implemented according to the present invention is capable of classifying and detecting the pain status in a patient as shown in Tables 6, 8 and 10

REFERENCES

(97) Akselrod, S., D. Gordon, F. A. Ubel, D. C. Shannon, A. C. Berger, and R. J. Cohen. “Power spectrum analysis of heart rate fluctuation: a quantitative probe of beat-to-beat cardiovascular control.” Science 213, no. 4504 (July 1981): 220-222.

(98) Belkin, M., and P. Niyogi. “Laplacian eigenmaps for dimensionality reduction and data representation.” Neural Computation 15 (2003): 13731396.

(99) Bishop, C. M. Pattern recognition and machine learning. New York: Springer, 2006.

(100) Borg, I., and P. Groenen. Modern Multidimensional Scaling: theory and applications. New York: Springer-Verlag, 2005.

(101) Breiman L, Friedman J, Olshen R, Stone C. Classification and Regression Trees. Belmont, Calif.: Wadsworth; 1984.

(102) Breiman, L, “Bagging Predictors.” Machine Learning 24, no. 2 (1996): 123-140.

(103) Breiman, L. “Random Forests.” Machine Learning V45, no. 1 (October 2001): 5-32.

(104) Breiman L. Random Forests, Machine Learning J., 2001b; 45:5-32.

(105) Bruhn, J., H. Ropcke, and A. Hoeft. “Approximate entropy as an electroencephalographic measure of anesthetic drug effect during desflurane anesthesia.” Anesthesiology 92, no. 3 (March 2000): 715-726.

(106) Chapelle, O., V. Vapnik, O. Bousquet, and S. Mukherjee. “Choosing Kernel Parameters for Support Vector Machines.” Machine Learning, 2000

(107) Cohen, D. H., and S. M. Sherman. “The autonomic nervous system and its central control.” Edited by R. Berne, and M. Levy, 322. St. Louis: C. V. Mosby, 1983.

(108) Coifman, R. R., et al. “Geometric diffusions as a tool for harmonic analysis and structure definition of data.” Proceedings of the National Academy of Sciences 102 (2005): 7426-7438.

(109) d'Aspremont, A., L. El Ghaoui, M. I. Jordan, and G. R. G. Lanckriet. “A direct formulation for sparse PCA using semidefinite programming.” Cambridge, Mass.: MIT Press, 2005.

(110) Deschamps, A., I. Kaufman, S. B. Backman, and G. Plourde. “Autonomic nervous system response to epidural analgesia in laboring patients by wavelet transform of heart rate and blood pressure variability.” Anesthesiology 101, no. 1 (July 2004): 21-27.

(111) Duda R, Hart P, Stork D. Pattern Classification. New York: Willey; 2000.

(112) Donoho, D. L., and C. Grimes. “Hessian eigenmaps: Locally linear embedding techniques for high-dimensional data.” PNAS 100, no. 10 (May 2003): 5591-5596.

(113) Goncharova, I. I., and J. S. Barlow. “Changes in EEG mean frequency and spectral purity during spontaneous alpha blocking.” Electroencephalography and clinical neurophysiology 76, no. 3 (September 1990): 197-204.

(114) Gronwall D, Wrightson P. Delayed Recovery of intellectual function after minor head injury. Lancet 2 (1974):605-609.

(115) Guignard, B. “Monitoring analgesia.” Best practice and research. Clinical anaesthesiology 20, no. 1 (March 2006): 161-180.

(116) Guo, Y., T. Hastie, and R. Tibshirani. “Regularized linear discriminant analysis and its application in microarrays.” Biostatistics (Oxford University Press) 8, no. 1 (January 2007): 86-100.

(117) Guyton, A. C. Human Physiology and Mechanisms of Disease. Philadelphia: Saunders, 1982.

(118) Fagrell, B. Microcirculation of the skin. The physiology and pharmacology of the microcirculation. Vol. 2. New York: Academic Press, 1984.

(119) Hastie, T., R. Tibshirani, and J. Friedman. The elements of statistical learning: data mining, inference, and prediction. New York: Springer, 2009.

(120) Hjorth, B. “The physical significance of time domain descriptors in EEG analysis.” Electroencephalography and clinical neurophysiology 34, no. 3 (March 1973): 321-325.

(121) Kohavi, R., and G. H. John. “Wrappers for Feature Subset Selection.” Artificial Intelligence 97, no. 1-2 (1997): 273324.

(122) Lidberg, L., and G. B. Wallin. “Sympathetic Skin Nerve Discharges in Relation to Amplitude of Skin Resistance Responses.” Psychophysiology 18, no. 3 (1981): 268-270.

(123) Malik, M. “Heart rate variability: standards of measurement, physiological interpretation and clinical use. Task Force of the European Society of Cardiology and the North American Society of Pacing and Electrophysiology.” Circulation 93, no. 5 (March 1996): 1043-1065.

(124) Moghaddam, B., Y. Weiss, and S. Avidan. “Spectral Bounds for Sparse PCA: Exact and Greedy Algorithms.” Edited by Weiss, Y., Scholkopf, B. and Platt, J., 915-922. Cambridge, Mass.: MIT Press, 2006.

(125) Moghaddam, B., Y. Weiss, and S. Avidan. “Generalized Spectral Bounds for Sparse LDA.” ICML. 2006a. 641-648.

(126) Oppenheim A V, Schafer R W. Discrete-time signal Processing. Upper Saddle River, N.J.: Prentice Hall; 1999.

(127) Pagani, M., O. Rimoldi, and A. Malliani. “Low-frequency components of cardiovascular variabilities as markers of sympathetic modulation.” Trends in pharmacological sciences 13, no. 2 (February 1992): 50-54.

(128) Pan, J., and W. J. Tompkins. “A real-time QRS detection algorithm.” IEEE Trans Biomed Eng 32, no. 3 (March 1985): 230-236.

(129) Pop-Jordanova, N., and J. Pop-Jordanov. “Spectrum-weighted EEG frequency (“brain-rate”) as a quantitative indicator of mental arousal.” Prilozi Makedonska akademija na naukite i umetnostite, Oddelenie za biolovski i medicinski nauki 26, no. 2 (December 2005): 35-42.

(130) Rangayyan R M. Biomedical signal analysis: A case study approach. IEEE Press, 2002.

(131) Sanei S, Chambers J A. EEG signal processing. New York: John Wiley & Sons; 2007.

(132) Schlogl, Alois. “A comparison of multivariate autoregressive estimators.” Signal Process. (Elsevier North-Holland, Inc.) 86, no. 9 (September 2006): 2426-2429

(133) Scholkopf, B., A. Smola, and K.-R. Muller. “Nonlinear component analysis as a kernel eigenvalue problem.” Neural Comput. (MIT Press) 10, no. 5 (1998): 1299-1319.

(134) Seitsonen, E. R., I. K. Korhonen, M. J. Van Gils, J. M. Lotjonen, K. T. Kortilla, and A. M. Yli-Hankala. “EEG spectral entropy, heart rate, photoplethysmography and motor responses to skin incision during sevoflurane anaesthesia.” Acta Anaesthesiol Scand, 49, (2005): 284292.

(135) Roweis, S., and L. Saul. “Nonlinear dimensionality reduction by locally linear embedding.” Science 290. (2000): 2323-2326.

(136) Tenenbaum, J., V. de Silva, and J. C. Langford. “A global geometric framework for nonlinear dimensionality reduction.” Science 290 (2000): 23192323.

(137) Thomas, D. W., and J. M. Evans. “Lower oesophageal contractility monitoring during anaesthesia for cardiac surgery: preliminary observations.” Annals of the Royal College of Surgeons of England 71, no. 5 (September 1989): 311-315.

(138) Tibshirani, R., T. Hastie, B. Narasimhan, and G. Chu. “Diagnosis of multiple cancer types by shrunken centroids of gene expression.” PNAS 99, no. 10 (May 2002): 6567-6572.

(139) van den Berg, A. A., D. Savva, and N. M. Honjol. “Attenuation of the haemodynamic responses to noxious stimuli in patients undergoing cataract surgery. A comparison of magnesium sulphate, esmolol, lignocaine, nitroglycerine and placebo given i.v. with induction of anaesthesia.” European Journal of Anaesthesiology 14, no. 02 (2006): 134-147.

(140) Vapnik, V. N. Statistical Learning Theory. Edited by Haykin S. New York: Wiley, 1998.

(141) Wackermann, J. “Towards a quantitative characterisation of functional states of the brain: from the non-linear methodology to the global linear description.” International Journal of Psychophysiology 34, no. 1 (October 1999): 65-80.

(142) Weiss, T., A. Del Bo, N. Reichek, and K. Engelman. “Pulse transit time in the analysis of autonomic nervous system effects on the cardiovascular system.” Psychophysiology 17, no. 2 (March 1980): 202-207.

(143) Zou, H., T. Hastie, and R. Tibshirani. “Sparse Principal Component Analysis.” Journal of Computational Graphical Statistics 15 (2006): 265286.

(144) While the invention has been described with respect to a limited number of embodiments, it will be appreciated that many variations, modifications and other applications of the invention may be made.