A Multimodal Deep Learning Approach to Intraoperative Nociception Monitoring: Integrating Electroencephalogram, Photoplethysmography, and Electrocardiogram
Pain and nociception are often confused, but they are distinct. Pain is a conscious, subjective experience . Under general anesthesia, however, “pain” is more accurately referred to as “nociception”, which reflects the physiological response to noxious stimuli in the absence of conscious perception. In the last two decades, significant progress has been made in developing methods to monitor nociception during anesthesia, particularly through algorithms that assess the nociception–antinociception (NAN) balance based on hemodynamic (e.g., ECG, PPG) and neural (e.g., EEG) signals. An increase in the sympathetic tone might be perceived as a reaction to painful stimuli , while conversely, an increase in the parasympathetic tone might be related to analgesia. The Surgical Pleth Index (SPI) is derived from PPG signals and used to monitor NAN; it uses a combination of normalized heartbeat intervals (HBI Norm ) and plethysmography amplitude (PPGA Norm ) . SPI is effective for nociception monitoring; however, since it is based on autonomic nervous system (ANS) activity, it might be confounded by other factors such as arrhythmia or antihypertensive drugs . The analgesia nociception index (ANI), another measure used to monitor NAN, is derived from high-frequency heartbeat intervals . The ANI is based on parasympathetic reactions and reflects noxious stimulations; however, it is also prone to arrhythmia and affected by parasympathetic reactions not related to pain. Moreover, the ANI was found to have a large inter-individual variability . In general anesthesia, EEG signals are used to predict the consciousness level; however, some studies show that EEG-derived indexes (e.g., BIS, SE, and RE) increased significantly in painful stimulus events , and thus, it might be helpful to provide information about nociception. qNOX is another index that is derived from EEG signals . It was developed using an Adaptive Neuro Fuzzy Inference System (ANFIS) by fitting different EEG frequency bands into a reference scale to predict noxious stimuli responses. Other methods have used single or a combination of parameters to monitor nociception . Pain/nociception is a subjective experience . Even with the absence of a clinical response, brain and spinal activities persist . Pain perception starts when nociceptors sense painful stimuli. Nociceptors send pain signals through the spinal cord to the brainstem, where the first stage of processing takes place. These signals then move upward to the Thalamus, which serves as a central relay station for sensory information. Further processing and relay occur in higher brain regions, including the frontal cortex, specifically the somatosensory cortex and anterior cingulate cortex . EEG bands were reported to react to several kinds of painful stimuli. It was found that EEG patterns are altered by noxious procedures in children, with local anesthesia mitigating this response . Additionally, under sevoflurane anesthesia, the neuromuscular block affects BIS and EEG responses to noxious electrical stimulation, suggesting genuine arousal differences . Clinical observations during tracheal intubation highlight alpha power changes as a sensitive indicator of cerebral activity modulation . Moreover, noxious stimulation affects cortical electrical activity levels, as shown by EEG parameters . The model selection process is challenging. Different models and methods were used in previous studies. For instance, SPI uses a linear model , while qNOX was developed using an ANFIS model , and the nociception level index was developed using random forest and linear regression models . Long Short-Term Memory Neural Network (LSTM) models  are known for their ability to retain important information over a long period; therefore, they have the advantage of using previous knowledge to make predictions efficiently. On the other hand, Multilayer Perceptron (MLP) models do not need previous knowledge , and yet they can predict nonlinear relationships and find hidden patterns within the data . While previous methods for nociception monitoring have demonstrated utility, the limitations of single-source approaches, particularly those reliant solely on the autonomic nervous system (ANS), prompt consideration of multimodal methods. ANS-based models, although informative, may be susceptible to non-nociceptive factors, potentially compromising the accuracy of nociception assessments. Moreover, the use of EEG data in isolation may not provide a comprehensive understanding of nociceptive states. In contrast, multimodal monitoring could offer the advantage of encompassing a broad spectrum of biomedical and clinical information. By integrating multiple physiological signals, such as EEG, PPG, and ECG, multimodal models have the capacity to capture diverse aspects of pain processing. Based on these considerations, we aim to address this challenge by implementing a multimodal approach utilizing EEG, PPG, and ECG signals to predict nociception during general anesthesia by integrating these three signals into one model. This study uses expert anesthesiologists’ nociception assessments (NOAs) based on clinical signs noted carefully during surgeries. 

2.1. Data Collection This work received approval from the institutional review board of the National Taiwan University Hospital (NTUH), and informed consent was obtained from all participants involved in this study. ECG, EEG, and PPG signals were collected using a Philips IntelliVue MP60 physiological signal monitor (Koninklijke Philips N.V, Amsterdam, The Netherlands), and the data were stored on a personal computer for further analysis. This study initially enrolled 142 patients from NTUH scheduled for general anesthesia surgery with inhalation anesthesia, during which propofol and fentanyl were administered. Another dataset from ECKH ( n = 10) was included in the analysis when it became available, the signals were collected using CARESCAPE B650 monitor (GE Healthcare Finland Oy, Helsinki, Finland); the anesthesia method was also inhalation anesthesia. Due to missing signal data in the NTUH dataset, the final analysis included 90 patients from NTUH and 10 patients from ECKH. The patients’ demographic data are listed in . Currently, no universally accepted gold standard for nociception exists, making assessment more challenging. The NOAs were scored based on each patient’s anesthesia record. The assessments were later digitized and stored in text files. After the surgeries were finished, each expert anesthesiologist (five in total) was appointed to assess the pain based on the information contained in the anesthesia record (e.g., HR, BP, anesthetic gas concentration, opioid dosages, etc.), and events, such as the eyelash reflex and laryngeal mask airway, were noted. As each anesthesiologist assessed pain individually, we synchronized the pain assessment data collected from multiple doctors and aligned them with the corresponding ECG, PPG, and EEG signals for analysis. Despite potential variations in assessment start times in a few cases, we matched data durations by comparing each patient’s assessment times and excluding offsets. The synchronization algorithm identified the latest start time (L) and the earliest end time (U) across all assessments, defining the analysis time interval [L, U]. Data points outside this interval were removed, ensuring only synchronized data were retained. Similar procedures were applied to align ECG, PPG, and EEG signals with synchronized pain score data. This process ensured consistency and the alignment of data across assessments. Our study used different parameters extracted from each signal. The ECG signals were used to extract the high-frequency waveform of heartbeats (RRHF) and its spectral power (RRHF PS ). The photoplethysmographic pulse wave amplitude (PPGA) and the area under the curve of PPG signals (PPG AUC ) were also extracted. Additionally, EEG signals were used to calculate the spectral power for each frequency band (i.e., delta, theta, alpha, beta, and gamma). The proposed model in this study is described in the flowchart shown in . All processing steps were conducted using MATLAB (version R2022b, MathWorks Inc., Natick, MA, USA) with the Signal Processing Toolbox. 2.2. Signal Processing and Feature Extraction First, RR series were extracted from ECG (512 Hz) using Pan–Tompkin’s algorithm . Then, we applied wavelet decomposition to the raw data and reconstructed only the RRHF in the range of 0.15–0.5 Hz, which is related to parasympathetic activity. Moreover, we calculated the power spectrum for the RRHF series over 64 s segments with a Kaiser window of 16 s. The PPG signals with a 128 Hz sampling frequency were filtered as needed. We used a Chebyshev II fifth-order filter with a cutoff of 0.5–8 Hz. For peak detection, we used the automatic multiscale-based peak detection algorithm . The long-term trend of PPGA was more important than the instant value changes; therefore, we applied a moving average filter of a 16 s window on the raw data. The PPG AUC was calculated for each 64 s window. shows a bandpass-filtered ECG window with its corresponding R peaks detected, and shows a PPG 64 s window before and after filtering. EEG (128 Hz) was bandpass filtered between 0.5 and 48 Hz to keep the signal of interest that covered all the EEG bands. Then, the power spectrum was calculated by thresholding over each band’s frequency limits. An example of the raw EEG and its 5 bands are shown in . The EEG signal was first divided into 64 s segments with a 5 s time step; then, for each segment, the power was calculated using a Kaiser window of 16 s with no overlapping. shows the procedure for data windowing. The resulting power array was then averaged to obtain one value for each 64 s window. Although this is not standard, the used limits were as follows: Delta: 0.5–3.5. Theta: 4–7.5. Alpha: 8–12. Beta: 13–30. Gamma: 30.5–48. 2.3. Normalization Initially, the histogram normalization described in  was used with a few adjustments. First, the standard deviation for the group dataset was not fixed; instead, it was the actual measure of the dataset. Second, the original distribution parameters were not combined; however, the transformation was calculated based on two distributions. The first transformation was based on the group dataset’s cumulative distribution function (cdf), and the second transformation was based on an accumulated individual dataset. At the beginning of the surgery, the weight for the group transformation was 100%. As more data became available, the weight for group data was reduced linearly until 10 min of the individual data was accumulated. The weight after that was fixed at 0.7 for individuals and 0.3 for groups. The transformations were calculated as: (1) X n o r m = w · F x x t + ( 1 − w ) · F y ( x t ) , where X norm is the normalized value of x i j . F x ( x ( t ) ) is the cumulative probability of x ( t ) in the individual dataset (from the beginning of the surgery until the current value). F y ( x ( t ) ) is the cumulative probability of x ( t ) in the group dataset. w is the weight for the individual dataset, adjusted based on t , which is the window number. w starts at 0 and smoothly adjusts to 0.7 over the segments. The histogram normalization is modeled in . When the data were normalized using histogram normalization, we noticed that the correlation between the normalized features and the pain assessment might be affected. Moreover, the variability in the original features was negatively affected in some cases; therefore, we decided to use Min–Max and z-score normalizations. The implementation was similar to histogram normalization, where the data were separated into group and individual datasets and the weights changed linearly based on the collected individual data from 100% to 30% for the group dataset. The z-score was previously utilized for online data normalization , and the equation in our study is as follows: (2) X n o r m = w 1 · x t − m e a n x a c c u m u l a t e d s t d x a c c u m u l a t e d + w 2 · x t − m e a n x g r o u p s t d x g r o u p , where X norm is the normalized value of x ( t ) . x a c c u m u l a t e d is the collected data from the beginning of the surgery to the current time. x g r o u p is the group dataset. For Min–Max normalization, we followed a similar procedure as implemented in ; the used equation is: (3) X n o r m = w 1 · x t − m i n x a c c u m u l a t e d m a x x a c c u m u l a t e d − m i n x a c c u m u l a t e d + w 2 · x t − m i n x g r o u p m a x x g r o u p − m i n x g r o u p , where X norm is the normalized value of x ( t ) . x a c c u m u l a t e d is the collected data from the beginning of the surgery to the current time. x g r o u p is the group dataset. To compare the three methods, the residuals and correlation between the original parameters and their transformed values were compared. As the data were available offline, we assumed that a perfect normalization would be achieved by normalizing the data while minimizing both intra- and inter-patient variability. Therefore, the normalization of the offline data to obtain a perfect normalization was fixed with a weight of 0.7 for the individual data, as we already had the data, and with a weight of 0.3 for the group data. The obtained values were then compared to the online implementation of the normalization. The comparison was performed by calculating the MAE and correlation between the normalized data. The best method was selected based on a metric score; the correlation and MAE were normalized and given a percentage based on their values. As both MAE and correlation were equally important, the score was calculated as follows: (4) S c o r e  = C o r r e l a t i o n n o r m a l i z e d + M A E n o r m a l i z e d 2 2.4. Deep Learning Training The models’ parameters were selected after applying an optimizable approach. The used models were MLP and LSTM. The models hyperparameters are listed in . The input parameters for the model were the five spectral power variables extracted from EEG signals, the spectral power of RRHF series, and the PPGA and PPG AUC series. We first assumed that the model did not need to be complex; however, to cover every possibility, we set the number of hidden layers to be between 2 and 10. The learning rate was selected between 0.1 and 0.0001, and the batch size was random between 64 and 512. The selected optimal parameters are described in . and show the models’ architecture. 2.5. Statistical Analysis Three tests were applied to measure the agreement between each doctor’s assessment. Analysis of variance (ANOVA) with a 0.05 significance level was used to test the similarity in the distributions and variance. To ensure consistency and agreement between the experts’ assessments, we applied the Bland–Altman test  and intraclass correlation (ICC)  with a 0.05 significance level. The ICC test was performed on individual assessments and on the whole data. The correlation coefficients and t -tests were calculated for pain assessments and the parameters used in this study. The used models were evaluated using the mean squared error (MSE) and R 2 on the training and validation datasets. The test (10 patients) dataset was evaluated using the mean absolute error (MAE) and correlation analysis. For the ECKH data of 10 patients, we compared the differences of the mean at three different surgical stages: Intubation (t 1 ), Incision (t 2 ), and Extubation (t 3 ). The receiver operating characteristic curve (ROC) was used to evaluate the models’ output in predicting nociception at the defined events for each patient individually (ECKH, n = 10). 

This work received approval from the institutional review board of the National Taiwan University Hospital (NTUH), and informed consent was obtained from all participants involved in this study. ECG, EEG, and PPG signals were collected using a Philips IntelliVue MP60 physiological signal monitor (Koninklijke Philips N.V, Amsterdam, The Netherlands), and the data were stored on a personal computer for further analysis. This study initially enrolled 142 patients from NTUH scheduled for general anesthesia surgery with inhalation anesthesia, during which propofol and fentanyl were administered. Another dataset from ECKH ( n = 10) was included in the analysis when it became available, the signals were collected using CARESCAPE B650 monitor (GE Healthcare Finland Oy, Helsinki, Finland); the anesthesia method was also inhalation anesthesia. Due to missing signal data in the NTUH dataset, the final analysis included 90 patients from NTUH and 10 patients from ECKH. The patients’ demographic data are listed in . Currently, no universally accepted gold standard for nociception exists, making assessment more challenging. The NOAs were scored based on each patient’s anesthesia record. The assessments were later digitized and stored in text files. After the surgeries were finished, each expert anesthesiologist (five in total) was appointed to assess the pain based on the information contained in the anesthesia record (e.g., HR, BP, anesthetic gas concentration, opioid dosages, etc.), and events, such as the eyelash reflex and laryngeal mask airway, were noted. As each anesthesiologist assessed pain individually, we synchronized the pain assessment data collected from multiple doctors and aligned them with the corresponding ECG, PPG, and EEG signals for analysis. Despite potential variations in assessment start times in a few cases, we matched data durations by comparing each patient’s assessment times and excluding offsets. The synchronization algorithm identified the latest start time (L) and the earliest end time (U) across all assessments, defining the analysis time interval [L, U]. Data points outside this interval were removed, ensuring only synchronized data were retained. Similar procedures were applied to align ECG, PPG, and EEG signals with synchronized pain score data. This process ensured consistency and the alignment of data across assessments. Our study used different parameters extracted from each signal. The ECG signals were used to extract the high-frequency waveform of heartbeats (RRHF) and its spectral power (RRHF PS ). The photoplethysmographic pulse wave amplitude (PPGA) and the area under the curve of PPG signals (PPG AUC ) were also extracted. Additionally, EEG signals were used to calculate the spectral power for each frequency band (i.e., delta, theta, alpha, beta, and gamma). The proposed model in this study is described in the flowchart shown in . All processing steps were conducted using MATLAB (version R2022b, MathWorks Inc., Natick, MA, USA) with the Signal Processing Toolbox. 

First, RR series were extracted from ECG (512 Hz) using Pan–Tompkin’s algorithm . Then, we applied wavelet decomposition to the raw data and reconstructed only the RRHF in the range of 0.15–0.5 Hz, which is related to parasympathetic activity. Moreover, we calculated the power spectrum for the RRHF series over 64 s segments with a Kaiser window of 16 s. The PPG signals with a 128 Hz sampling frequency were filtered as needed. We used a Chebyshev II fifth-order filter with a cutoff of 0.5–8 Hz. For peak detection, we used the automatic multiscale-based peak detection algorithm . The long-term trend of PPGA was more important than the instant value changes; therefore, we applied a moving average filter of a 16 s window on the raw data. The PPG AUC was calculated for each 64 s window. shows a bandpass-filtered ECG window with its corresponding R peaks detected, and shows a PPG 64 s window before and after filtering. EEG (128 Hz) was bandpass filtered between 0.5 and 48 Hz to keep the signal of interest that covered all the EEG bands. Then, the power spectrum was calculated by thresholding over each band’s frequency limits. An example of the raw EEG and its 5 bands are shown in . The EEG signal was first divided into 64 s segments with a 5 s time step; then, for each segment, the power was calculated using a Kaiser window of 16 s with no overlapping. shows the procedure for data windowing. The resulting power array was then averaged to obtain one value for each 64 s window. Although this is not standard, the used limits were as follows: Delta: 0.5–3.5. Theta: 4–7.5. Alpha: 8–12. Beta: 13–30. Gamma: 30.5–48. 

Initially, the histogram normalization described in  was used with a few adjustments. First, the standard deviation for the group dataset was not fixed; instead, it was the actual measure of the dataset. Second, the original distribution parameters were not combined; however, the transformation was calculated based on two distributions. The first transformation was based on the group dataset’s cumulative distribution function (cdf), and the second transformation was based on an accumulated individual dataset. At the beginning of the surgery, the weight for the group transformation was 100%. As more data became available, the weight for group data was reduced linearly until 10 min of the individual data was accumulated. The weight after that was fixed at 0.7 for individuals and 0.3 for groups. The transformations were calculated as: (1) X n o r m = w · F x x t + ( 1 − w ) · F y ( x t ) , where X norm is the normalized value of x i j . F x ( x ( t ) ) is the cumulative probability of x ( t ) in the individual dataset (from the beginning of the surgery until the current value). F y ( x ( t ) ) is the cumulative probability of x ( t ) in the group dataset. w is the weight for the individual dataset, adjusted based on t , which is the window number. w starts at 0 and smoothly adjusts to 0.7 over the segments. The histogram normalization is modeled in . When the data were normalized using histogram normalization, we noticed that the correlation between the normalized features and the pain assessment might be affected. Moreover, the variability in the original features was negatively affected in some cases; therefore, we decided to use Min–Max and z-score normalizations. The implementation was similar to histogram normalization, where the data were separated into group and individual datasets and the weights changed linearly based on the collected individual data from 100% to 30% for the group dataset. The z-score was previously utilized for online data normalization , and the equation in our study is as follows: (2) X n o r m = w 1 · x t − m e a n x a c c u m u l a t e d s t d x a c c u m u l a t e d + w 2 · x t − m e a n x g r o u p s t d x g r o u p , where X norm is the normalized value of x ( t ) . x a c c u m u l a t e d is the collected data from the beginning of the surgery to the current time. x g r o u p is the group dataset. For Min–Max normalization, we followed a similar procedure as implemented in ; the used equation is: (3) X n o r m = w 1 · x t − m i n x a c c u m u l a t e d m a x x a c c u m u l a t e d − m i n x a c c u m u l a t e d + w 2 · x t − m i n x g r o u p m a x x g r o u p − m i n x g r o u p , where X norm is the normalized value of x ( t ) . x a c c u m u l a t e d is the collected data from the beginning of the surgery to the current time. x g r o u p is the group dataset. To compare the three methods, the residuals and correlation between the original parameters and their transformed values were compared. As the data were available offline, we assumed that a perfect normalization would be achieved by normalizing the data while minimizing both intra- and inter-patient variability. Therefore, the normalization of the offline data to obtain a perfect normalization was fixed with a weight of 0.7 for the individual data, as we already had the data, and with a weight of 0.3 for the group data. The obtained values were then compared to the online implementation of the normalization. The comparison was performed by calculating the MAE and correlation between the normalized data. The best method was selected based on a metric score; the correlation and MAE were normalized and given a percentage based on their values. As both MAE and correlation were equally important, the score was calculated as follows: (4) S c o r e  = C o r r e l a t i o n n o r m a l i z e d + M A E n o r m a l i z e d 2 

The models’ parameters were selected after applying an optimizable approach. The used models were MLP and LSTM. The models hyperparameters are listed in . The input parameters for the model were the five spectral power variables extracted from EEG signals, the spectral power of RRHF series, and the PPGA and PPG AUC series. We first assumed that the model did not need to be complex; however, to cover every possibility, we set the number of hidden layers to be between 2 and 10. The learning rate was selected between 0.1 and 0.0001, and the batch size was random between 64 and 512. The selected optimal parameters are described in . and show the models’ architecture. 

Three tests were applied to measure the agreement between each doctor’s assessment. Analysis of variance (ANOVA) with a 0.05 significance level was used to test the similarity in the distributions and variance. To ensure consistency and agreement between the experts’ assessments, we applied the Bland–Altman test  and intraclass correlation (ICC)  with a 0.05 significance level. The ICC test was performed on individual assessments and on the whole data. The correlation coefficients and t -tests were calculated for pain assessments and the parameters used in this study. The used models were evaluated using the mean squared error (MSE) and R 2 on the training and validation datasets. The test (10 patients) dataset was evaluated using the mean absolute error (MAE) and correlation analysis. For the ECKH data of 10 patients, we compared the differences of the mean at three different surgical stages: Intubation (t 1 ), Incision (t 2 ), and Extubation (t 3 ). The receiver operating characteristic curve (ROC) was used to evaluate the models’ output in predicting nociception at the defined events for each patient individually (ECKH, n = 10). 

3.1. NOA The pain assessments from each doctor were processed carefully. As five anesthesiologists provided the assessments, we first used the ANOVA test. Doctor D provided the lowest assessment, and it was significantly different than the other four doctors . Therefore, we only used assessments from Doctors A, B, C, and E. Moreover, we applied Bland–Altman and ICC  tests. When the pain assessments from all four doctors were pulled together, the ICC was 0.73 with a confidence interval of [0.69 0.75]. The ICC was also applied to each patient’s assessments individually. shows a histogram of the ICC values among the patients, and shows the ANOVA test results. The agreement between the doctors’ assessments was evaluated using Bland–Altman analysis and pairwise comparisons, with the results shown in and . The analysis results were consistent across the doctors in terms of agreement, where all data fell within the 95% limits of agreement (LoA). The lowest bias (−0.7) was observed between Doctors B and C, indicating minimal differences in their assessments. The highest difference was found between Doctors C and E (bias = 6.25) and Doctors B and E (bias = 5.55). The standard deviation of the differences ranged from 7.4 to 9.1. 3.2. Parameter Correlations After we matched the NOA for the four doctors, the average was taken for further analysis. The Pearson correlation coefficient between the parameters and the average assessment were calculated for each patient individually. The results were pulled together and are shown in . PPG AUC had the highest individual positive correlation (0.85) and mean (0.46). The PPGA correlation had a mean value of −0.32; it was mostly negative, with the highest value of −0.82. The power calculations for EEG bands correlated negatively in most cases except for gamma, which had a close mean value (0.40) to PPG AUC . RRHF PS had the lowest correlation along with the beta band (−0.003 and −0.004 mean values, respectively). shows an example of the used parameters, where the curves show a good positive or negative correlation with the pain assessment. The correlation analysis for the same subject is shown in . 3.3. Parameter Normalization The normalization methods resulted in different MAE and correlation values. The lowest MAE was found for Min–Max normalization, and the highest correlation was found for histogram normalization. The results are shown in . An example of normalization is shown in . For this case, the lowest MAE (0.08) was associated with Min–Max normalization, and the highest correlation was 0.95 for z-score normalization. 3.4. Pain/Nociception Models Initially, we had 90 patients’ datasets. When more datasets from another hospital (ECKH-2023 datasets) became available, we decided to assess the model from a medical and engineering perspective. For NTUH (2015), the pain assessments were based on anesthetic records. However, these records were not easily available to us at the moment of this study; therefore, we were unable to assess the 10 patient test datasets from a medical perspective. Thus, the NTUH data were evaluated using MSE and R 2 on the training and validation datasets and correlation analysis and MAE on the test dataset. For ECKH, the needed information, such as painful stimuli events, was available; therefore, medical analysis was possible. 3.4.1. MLP Model The MLP losses were 0.008 and 0.007 on training and validation, respectively. The R 2 was 0.82 on training and 0.80 on validation. For the test dataset (NTUH), the mean values for correlation and MAE were 0.78 and 5.1, respectively. shows the curves for training and validation, and shows the distribution of the correlation and MAE values for the test datasets. For the 10 patients’ datasets from ECKH, the mean value before and after the three defined events was calculated. The mean values after t 1 increased from 48 to 57. While at t 2 , the mean value increased the least, from 41 to 50, and at t 3 , the mean value increased the most, from 48 to 62. shows the changes in the mean values before and after each event (upper plot) and the difference in the mean NOA at each event (lower plot). The upper plot shows how the mean values before and after each event varied across the 10 patients. The lower plot demonstrates the differences in the mean NOA at each event between the patients. The highest difference was found at t 3 , with a value of 13, and the lowest was at t 2 , with a value of 8. The difference in the NOA means at t 1 was also clearly noticed, with a value of 9. An example of the predicted pain/nociception from the two models is shown in (NTUH dataset). 3.4.2. LSTM Model The training and validation losses for the LSTM model were 0.012 and 0.011, respectively. The R 2 was 0.87 and 0.81 on the training and validation, respectively. The mean value of the correlation and MAE were 0.85 and 3.6, respectively. shows the trends of the training and validation. shows the correlation and MAE values for each patient. The results were consistent, except for patient 3, who had a high MAE of 7.8. For the 10 patients’ datasets from ECKH, the mean values before and after the three defined events increased less than the MLP model. The mean values before and after t 1 were 53 and 55, respectively, and it was the lowest. The mean values before and after t 2 were 47 and 51, respectively, and at t 3 , the values were 51 and 55 before and after, respectively. shows the changes in the mean values before and after each event (upper plot) and the difference in the mean NOA at each event (lower plot). The upper plot shows how the mean values before and after each event varied across the 10 patients. The lower plot demonstrates the differences in the mean NOA at each event between the patients. The highest difference was found at t 2 and t 3 , with a value of 4, and the lowest was at t 1 , with a value of 2.5. shows an example of the models’ predictions on the ECKH dataset. The ROC curve results are shown in . The results were obtained by classifying NOAs based on the events. The figure shows that the MLP model outperformed the LSTM model in predicting noxious stimuli. The mean value for MLP AUC was 0.90; the highest was 0.95, and the lowest was 0.82. The LSTM demonstrated low to moderate reactivity; the mean AUC value was 0.68, whereas the lowest and highest values were 0.44 and 0.82, respectively. 

The pain assessments from each doctor were processed carefully. As five anesthesiologists provided the assessments, we first used the ANOVA test. Doctor D provided the lowest assessment, and it was significantly different than the other four doctors . Therefore, we only used assessments from Doctors A, B, C, and E. Moreover, we applied Bland–Altman and ICC  tests. When the pain assessments from all four doctors were pulled together, the ICC was 0.73 with a confidence interval of [0.69 0.75]. The ICC was also applied to each patient’s assessments individually. shows a histogram of the ICC values among the patients, and shows the ANOVA test results. The agreement between the doctors’ assessments was evaluated using Bland–Altman analysis and pairwise comparisons, with the results shown in and . The analysis results were consistent across the doctors in terms of agreement, where all data fell within the 95% limits of agreement (LoA). The lowest bias (−0.7) was observed between Doctors B and C, indicating minimal differences in their assessments. The highest difference was found between Doctors C and E (bias = 6.25) and Doctors B and E (bias = 5.55). The standard deviation of the differences ranged from 7.4 to 9.1. 

After we matched the NOA for the four doctors, the average was taken for further analysis. The Pearson correlation coefficient between the parameters and the average assessment were calculated for each patient individually. The results were pulled together and are shown in . PPG AUC had the highest individual positive correlation (0.85) and mean (0.46). The PPGA correlation had a mean value of −0.32; it was mostly negative, with the highest value of −0.82. The power calculations for EEG bands correlated negatively in most cases except for gamma, which had a close mean value (0.40) to PPG AUC . RRHF PS had the lowest correlation along with the beta band (−0.003 and −0.004 mean values, respectively). shows an example of the used parameters, where the curves show a good positive or negative correlation with the pain assessment. The correlation analysis for the same subject is shown in . 

The normalization methods resulted in different MAE and correlation values. The lowest MAE was found for Min–Max normalization, and the highest correlation was found for histogram normalization. The results are shown in . An example of normalization is shown in . For this case, the lowest MAE (0.08) was associated with Min–Max normalization, and the highest correlation was 0.95 for z-score normalization. 

Initially, we had 90 patients’ datasets. When more datasets from another hospital (ECKH-2023 datasets) became available, we decided to assess the model from a medical and engineering perspective. For NTUH (2015), the pain assessments were based on anesthetic records. However, these records were not easily available to us at the moment of this study; therefore, we were unable to assess the 10 patient test datasets from a medical perspective. Thus, the NTUH data were evaluated using MSE and R 2 on the training and validation datasets and correlation analysis and MAE on the test dataset. For ECKH, the needed information, such as painful stimuli events, was available; therefore, medical analysis was possible. 3.4.1. MLP Model The MLP losses were 0.008 and 0.007 on training and validation, respectively. The R 2 was 0.82 on training and 0.80 on validation. For the test dataset (NTUH), the mean values for correlation and MAE were 0.78 and 5.1, respectively. shows the curves for training and validation, and shows the distribution of the correlation and MAE values for the test datasets. For the 10 patients’ datasets from ECKH, the mean value before and after the three defined events was calculated. The mean values after t 1 increased from 48 to 57. While at t 2 , the mean value increased the least, from 41 to 50, and at t 3 , the mean value increased the most, from 48 to 62. shows the changes in the mean values before and after each event (upper plot) and the difference in the mean NOA at each event (lower plot). The upper plot shows how the mean values before and after each event varied across the 10 patients. The lower plot demonstrates the differences in the mean NOA at each event between the patients. The highest difference was found at t 3 , with a value of 13, and the lowest was at t 2 , with a value of 8. The difference in the NOA means at t 1 was also clearly noticed, with a value of 9. An example of the predicted pain/nociception from the two models is shown in (NTUH dataset). 3.4.2. LSTM Model The training and validation losses for the LSTM model were 0.012 and 0.011, respectively. The R 2 was 0.87 and 0.81 on the training and validation, respectively. The mean value of the correlation and MAE were 0.85 and 3.6, respectively. shows the trends of the training and validation. shows the correlation and MAE values for each patient. The results were consistent, except for patient 3, who had a high MAE of 7.8. For the 10 patients’ datasets from ECKH, the mean values before and after the three defined events increased less than the MLP model. The mean values before and after t 1 were 53 and 55, respectively, and it was the lowest. The mean values before and after t 2 were 47 and 51, respectively, and at t 3 , the values were 51 and 55 before and after, respectively. shows the changes in the mean values before and after each event (upper plot) and the difference in the mean NOA at each event (lower plot). The upper plot shows how the mean values before and after each event varied across the 10 patients. The lower plot demonstrates the differences in the mean NOA at each event between the patients. The highest difference was found at t 2 and t 3 , with a value of 4, and the lowest was at t 1 , with a value of 2.5. shows an example of the models’ predictions on the ECKH dataset. The ROC curve results are shown in . The results were obtained by classifying NOAs based on the events. The figure shows that the MLP model outperformed the LSTM model in predicting noxious stimuli. The mean value for MLP AUC was 0.90; the highest was 0.95, and the lowest was 0.82. The LSTM demonstrated low to moderate reactivity; the mean AUC value was 0.68, whereas the lowest and highest values were 0.44 and 0.82, respectively. 

The MLP losses were 0.008 and 0.007 on training and validation, respectively. The R 2 was 0.82 on training and 0.80 on validation. For the test dataset (NTUH), the mean values for correlation and MAE were 0.78 and 5.1, respectively. shows the curves for training and validation, and shows the distribution of the correlation and MAE values for the test datasets. For the 10 patients’ datasets from ECKH, the mean value before and after the three defined events was calculated. The mean values after t 1 increased from 48 to 57. While at t 2 , the mean value increased the least, from 41 to 50, and at t 3 , the mean value increased the most, from 48 to 62. shows the changes in the mean values before and after each event (upper plot) and the difference in the mean NOA at each event (lower plot). The upper plot shows how the mean values before and after each event varied across the 10 patients. The lower plot demonstrates the differences in the mean NOA at each event between the patients. The highest difference was found at t 3 , with a value of 13, and the lowest was at t 2 , with a value of 8. The difference in the NOA means at t 1 was also clearly noticed, with a value of 9. An example of the predicted pain/nociception from the two models is shown in (NTUH dataset). 

The training and validation losses for the LSTM model were 0.012 and 0.011, respectively. The R 2 was 0.87 and 0.81 on the training and validation, respectively. The mean value of the correlation and MAE were 0.85 and 3.6, respectively. shows the trends of the training and validation. shows the correlation and MAE values for each patient. The results were consistent, except for patient 3, who had a high MAE of 7.8. For the 10 patients’ datasets from ECKH, the mean values before and after the three defined events increased less than the MLP model. The mean values before and after t 1 were 53 and 55, respectively, and it was the lowest. The mean values before and after t 2 were 47 and 51, respectively, and at t 3 , the values were 51 and 55 before and after, respectively. shows the changes in the mean values before and after each event (upper plot) and the difference in the mean NOA at each event (lower plot). The upper plot shows how the mean values before and after each event varied across the 10 patients. The lower plot demonstrates the differences in the mean NOA at each event between the patients. The highest difference was found at t 2 and t 3 , with a value of 4, and the lowest was at t 1 , with a value of 2.5. shows an example of the models’ predictions on the ECKH dataset. The ROC curve results are shown in . The results were obtained by classifying NOAs based on the events. The figure shows that the MLP model outperformed the LSTM model in predicting noxious stimuli. The mean value for MLP AUC was 0.90; the highest was 0.95, and the lowest was 0.82. The LSTM demonstrated low to moderate reactivity; the mean AUC value was 0.68, whereas the lowest and highest values were 0.44 and 0.82, respectively. 

Many indexes were developed specifically to monitor the nociception–antinociception balance during general anesthesia , and all were based on different parameters or a combination of those parameters (e.g., PPGA, HBI, EEG features). All of them have proven utility and were able to react to painful stimuli events. Single-source indexes were reported to react to other factors that do not relate to nociception . On the other hand, we believe that single-source indexes might ignore other information that could be helpful to assess NAN. Therefore, in our study, we developed two models to predict nociception in various stages of surgery. First, the ANOVA test was used to check the distributions and variance in the assessments among the doctors. Doctor D showed disagreement with the others; therefore, the related assessments were excluded. To further prove that the other doctors’ assessments were consistent and reliable, we applied ICC and Bland–Altman tests on the assessments. The ICC was applied to the data when pulled together and individually. In both cases, the results proved the reliability between the assessment, as shown in (0.73 on the whole data). The Bland–Altman test also showed that 95% of the data was in the upper and lower limits of agreement ( and ); therefore, the assessments were used in further analysis. The parameter selection was based on the clinical parameters that were previously used to assess NAN . PPGA and HBI were used for SPI development in a single-source model  and in a multimodal approach , while EEG parameters were used in other methods . A correlation analysis was applied to the selected parameters with the pain assessments, and the results proved the validity of our pain assessments. The SPI was calculated using the equation (SPI = 100 − (0.7 × PPGA + 0.3 × HBI), where PPGA has the most contribution. This equation implies a negative correlation between nociception and PPGA or HBI, which was the case in our analysis, where PPGA correlated negatively with the assessments and the values were distributed in the range of −(0.2–0.82), as shown in . This agrees with the correlation analysis performed in  with pain stimulus intensity. The EEG parameters also showed an agreement with the results in , where the alpha band was reported to drop in the presence of painful stimuli. Other studies showed that the delta  and beta bands  increased with painful stimuli, and in our case, this was noted in many cases. Inter- and intra-patient variability needs to be addressed carefully. Previously, methods such as histogram normalization and z-score normalization were adapted for online usage. In our study, both methods were effective; however, we noticed that some parameters’ properties were affected after applying those techniques. Therefore, we adapted another normalization method (i.e., Min–Max normalization) for our use. The comparison between the three methods showed that for our data, the Min–Max method is the most effective in retaining the original parameter properties while reducing both inter- and intra-patient variability. In our study, both models used information from many sources, and both were based on experts’ assessments of pain. The models were evaluated in two groups: the first group was used to evaluate the models’ performances with the pain assessments (NTUH), and the second group was used to evaluate the models’ performances in the operation room (ECKH). On the NTUH datasets, the MLP model had a higher MAE and lower correlation compared to the LSTM model ( and ). When the predictions were visualized, however, we noticed that both models had predictions that generalized well around the pain assessments, with values from the LSTM model being lower and smoother. This difference was clear, and the LSTM predictions were not conclusive whether the current predicted value was related to nociception or not. On the other hand, the MLP predictions were more precise about the current state, whether it was related to nociception or not. This was proven in the prediction from the ECKH datasets. The predicted values before and after surgical events were compared. For the MLP model , the mean value after t 1 and t 3 increased the most, and it was more spread out after t3 (52–69). The mean value of the difference before and after each event provided evidence that the MLP model was able to predict nociception at painful stimuli events, where it was 9 and 12 at t 1 and t 3 , respectively. While smaller than t 1 and t 3 , the increase after t 2 was still noticed; this was also proven in the difference before and after t 2 , where the mean value was 8. On the other hand, for the LSTM model , the mean values of the predicted nociception before and after each event were not large enough to distinguish between different levels of nociception. The largest increase in the predicted value was noticed after t2 and t3, where the difference between the means before and after the two events was also small (4). The obtained results in and provide insights into the models’ behaviors at the defined events. In particular, the difference in the NOA mean at each event describes the ability of the model to distinguish a noxious stimulus. This was further proven by applying ROC on the defined events . The MLP model demonstrated more accurate predictions, achieving a mean AUC of 0.90 (range 0.82–0.95). Hence, this indicates a distinctive assessment of nociception at various stages of surgery. In contrast, the LSTM model provided more variable performance, with a mean AUC value of 0.68 (range 0.44–0.82). Considering the metrics used in our study, there is a clear gap in the performance of the two models that was noticed when evaluating the ECKH dataset. This could be related to different factors, including the complexity of LSTM models. Our iterative approach in training the LSTM model did not improve the results. The selection for LSTM was based on the fact that it considers previous knowledge, which could be beneficial in some cases. In our case, however, this might have been one of the limitations to our model. Moreover, we hypothesize that the significant difference in the models’ performances could be related to the fact that MLP models could capture hidden patterns and complex relationships accurately. Delta arousal (power increase) and alpha dropout (power decrease) are related to noxious stimuli. In , we show how every single parameter reacted during surgery against MLP and LSTM predictions. At t 1 , MLP predictions increased, which corresponded to delta arousal and alpha dropout. At t 2 , less reactive behavior was noticed from alpha and delta; similarly, PPGA smoothly decreased, and changes were not noticed in PPG AUC . Although MLP predictions reacted to these features, resulting in higher predictions, the LSTM model was unable to react to these changes. At t 3 , small changes in delta and alpha and a slow decrease in PPGA were noticed. Unlike the LSTM model, the MLP model was able to react to these changes. These results highlight the need for multiple parameter approaches. A model utilizing one parameter might ignore other information. An increase in PPGA is related to analgesia. As seen in the figure between t 1 and t 2 , PPGA increased in this time segment; this increase could be related to other factors. In the absence of other parameters, a model would predict this value as a decrease in nociception; however, in our case, the usage of multiple parameters provided sufficient data to the MLP model, which identified this as a noxious stimulus. Our study aimed to enhance nociception monitoring in general anesthesia by developing a multimodal framework that integrates multiple physiological parameters, offering more comprehensive information for nociception assessment. We applied a comparative analysis of two models (MLP and LSTM) for nociception prediction, demonstrating their relevant strengths and limitations in different clinical settings. Our evaluation across two hospitals further supports the generalizability of these findings and the ability of our MLP model to be used in different clinical settings and patient populations. Our analysis for the used signals demonstrated the capability of EEG signals as an important source for pain assessment, which was further validated with the models’ results. Additionally, utilizing Min–Max normalization for online and offline usage provided better results than histogram and z-score normalization. This was found when the correlation analysis was applied to the pain assessments with the used parameters. Initially, the raw parameters correlated to pain assessments; however, histogram normalization affected these results negatively. Exploring other methods used in the literature, such as the z-score, led to using another common normalization method (i.e., Min–Max normalization). Our analysis for model development focused on traditional regression metrics (MSE, MAE, and R 2 ), given the continuous nature of our datasets. For model validations, a statistical comparison in the mean difference between events was utilized, and the models’ abilities to differentiate nociceptive events were tested using ROC curves. However, we acknowledge the importance of including performance metrics such as F1, precision, and recall. These metrics could be adopted in the future with larger datasets and various anesthesia settings. Our previous studies about pain/nociception did not achieve this performance. We used SPI  and ANI  features to predict pain scores using LSTM and MLP models. The best model  had a correlation of 0.63 and an MAE of 3.42; however, despite this difference in performance, the predictions did not generalize well to the pain assessments, which was due to insufficient information provided to the model. Hence, we solved this in our present study, and the results improved. The types of surgeries and anesthesia administration methods could have an effect on the way the model predicts nociception. Our datasets were a combination of inhaled anesthesia and, in a few cases, TCI anesthesia. In both cases, propofol was used intravenously. Therefore, the findings in this study are more related to inhaled anesthesia. The MLP model was more useful in distinguishing the nociceptive reactions at different surgical events. However, this was noticed only when the models were evaluated on ECKH datasets. Considering the size of the dataset (10 patients), the size of the test dataset might not be enough to prove this. Therefore, future studies are needed to validate this result. 

We developed two multimodal models for nociception prediction in general anesthesia. The models use physiological parameters from various sources rather than relying solely on one signal’s information. The developed models were correlated with an expert assessment of pain that was based on clinical information related to nociception and surgical stimulation events. Our investigation of normalization techniques revealed that Min–Max normalization was optimal for both online and offline applications, preserving essential parameter characteristics. The models were evaluated across two hospitals with varying results. In the first group, both models demonstrated effective nociception prediction, with the LSTM model providing smoother predictions. In the second group, while LSTM models are typically powerful for sequence data problems, this advantage did not translate to better performance in our specific application. The MLP model successfully differentiated between nociception levels at various events, while the LSTM model’s predictions were less definitive. Our results are promising; however, further studies are needed to validate these findings in different anesthesia settings and a larger patient population. 

