The investigation was granted approval by the local Ethics Committee of Istituto Nazionale Tumori, Fondazione Pascale, Naples, Italy (protocol code 41/20 Oss; approval date: 26 November 2020), and all participants provided written informed consent. All methods were performed following the ethical standards as laid down in the Declaration of Helsinki and its later amendments or comparable ethical standards and it was registered with ClinicalTrials.gov, number (NCT04726228). The study employed a prospective observational design, utilizing medical records as its primary data source. The cohort for this study comprised adult patients undergoing treatment for cancer-related pain at the Istituto Nazionale Tumori, Fondazione Pascale, Italy. Clinical-instrumental assessments and collection of anamnestic data were conducted during a single outpatient visit at the pain therapy clinic. The process of recording biosignals was executed spanning a timeframe of approximately 5 min.

The analysis involved the examination of numerous variables, including demographic factors such as age, and body mass index (BMI). Clinical measurements included the number of comorbidities (categorized as “None” and at least one), Eastern Cooperative Oncology Group (ECOG) performance status (categorized as low ≤ 2, high > 2), metastases (yes/no), bone metastases (yes/no). Other clinical data concerned the type of anticancer therapy (chemotherapy, immunotherapy, surgery, and radiotherapy). Pain-related variables were pain type (nociceptive, neuropathic) and intensity (0–10 NRS), breakthrough cancer pain (BTcP) [19], and its features (type: nociceptive or neuropathic; 0–10 NRS intensity). Pain therapy included opioids, calculated as morphine equivalent dose (MED) for background pain, rapid onset opioids for BTcP, and pain adjuvants (corticosteroids, antidepressants, anticonvulsants, muscle relaxants, and anti-anxiety medications).

The objective variables were obtained from the analysis of the EDA and ECG biosignals. EDA reflects the sympathetic nervous system activity and serves as a valuable indicator for evaluating pain-induced neurocognitive stress. It detects alterations in the skin's electrical properties, triggered by sweat gland activation, ultimately resulting in increased skin conductance [20]. The continuous variations in skin conductance are termed the Skin Conductance Level (SCL), while the rapid responses occurring within seconds are identified as the Galvanic Skin Response (GSR). Both SCL and GSR play roles in the tonic and phasic components (Fig. 1).

Phasic and tonic components of the electrodermal activity (EDA) signal over time. Legend: The signal is recorded in microsiemens (µS)

For acquiring and analyzing biosignals, we employed a previously validated approach [21]. Regarding EDA biosignals, two EDA-derived analyses have been performed including the continuous decomposition analysis (CDA), and the trough-to-peak (TTP) analysis, both conducted in accordance to [22,23,24]. For the analysis, the mean amplitude of all above-threshold GSRs was calculated for both CDA and TTP analyses and adopted as a representative parameter of the EDA signal since it provides an overall synthetic measure of the average amplitude of registered skin conductance responses.

Concerning the ECG signal, temporal changes in inter-beat intervals offer a gauge of heart rate variability (HRV), closely associated with the autonomic nervous system (ANS) activity. These fluctuations can also signal physiological reactions to stressful or painful situations. Following our previous study [21], we focused on the R-R series of interbeat intervals, defined as the time between successive R waves of the QRS complex on the ECG waveform, [25], and then used it to extract time-domain parameters of the HRV. In particular, as reported in [21], R peaks were detected using a modified version of the Pan-Tompkins technique based on the ECG envelope and subsequent flattening of the signal to enhance the QRS-complexes as proposed in [26,27,28], and the raw series of differences between consecutive peaks was obtained (R-R). Then, using a recursive filtering process, the ECG-derived R-R time series were filtered to remove outliers and intervals that deviated most from the mean of the nearby RR intervals. Subsequently, from the ECG-derived R-R time series (Fig. 2), the heart rate (HR) time series has been obtained, and the standard deviation (SD) of both R-R and HR series has been calculated and used as a reference time-domain indicator of the variability of the heart rhythm.

Plot of the ECG and R-R intervals over time

As we described in a previous study [21], EDA and ECG signals have been acquired by using a BITalino device equipped with sensors for the recording of ECG and EDA signals. The BITalino platform is a hard-ware-affordable and open-source biosignals platform developed for physiological computing. Signals have been collected at a 1000 Hz sampling rate.

The sensitivity analysis tests were conducted using the Li and Yu [29] and Raikov test [30] on the continuous variables for assessing the potential Missing Not At Random (MNAR) process presence. Therefore, based on the information obtainable from the data and the context, the best conclusion on the mechanism generating the missing data will be reached. A Multiple Imputation analysis was finally performed on data to preserve the sample size.

A Multiple Factorial Analysis (MFA) was performed on data to detect multi-correlations and associations between numerical and categorical variables within the context of subjective and objective pain features. It is a statistical technique used for analyzing databases with large sets of variables that are suspected of being statistically related and describe the same set of observations or individuals [31]. This strategy is often implemented for addressing different types of variables (e.g., numerical, categorical) and for analyzing their relationships in a single analysis by creating principal components (PCs) which are assigned as a part of data behavior and variability. MFA allows analysts to explore patterns, similarities, and differences between the different sets of variables and observations [32]. In our analysis, a focus was applied to the main interest variables namely EDA and ECG derivates and BTcP presence (groups: No BTcP, BTcP), by reading 2D-projected modalities in the factorial plane respect with the BTcP intensity groups.

Within the realm of MFA, when dealing with complex multidimensional datasets featuring numerous categorical variables, the utilization of Multiple Correspondence Analysis (MCA) proves to be a valuable strategy. It enables the condensation of dataset dimensions into a selected set of categorical variables, effectively capturing the most essential information. Therefore, we adopted a data-driven approach known as Hierarchical Clustering on Principal Components (HCPC). This approach harmoniously integrates three established methods including MCA, or PC analysis (PCA) for numerical variables, hierarchical clustering, and the k-means algorithm, to derive a refined and enhanced cluster solution [33]. Through this amalgamation of techniques, we aimed to extract a more robust and accurate representation of patterns and relationships present within the data. Operatively, the HCPC method provides a PCA and, subsequently, a hierarchical clustering by implementing an agglomerative hierarchical tree; the ideal number of groups is obtained by a pruning approach; then k-means was performed on data by setting the number of groups, obtained as described above. Thus, different variables were inserted in the routine to find the BTcP-related risk groups. They included BTcP type, ECOG, metastases (No, Yes, bone metastases), chemotherapy, radiotherapy, and surgery, as well as drug therapy including MED (≤ 60 mg; > 60 mg) and adjuvants. Finally, a multivariable linear regression analysis was conducted to assess the main associations between MED (treated as a number and binary variable) and potentially related variables.

The data were analyzed using the R software version 4.2.3 (R Core Teams, R Foundation for Statistical Computing, Vienna, Austria). The toolkit included car, purr, boot, snow, misty, and naniar. The Mice package was adopted for the imputation of the missing data. FactoMineR was used for the implementation of factor analysis methods. Moreover, graphical packages were adopted for the visualization of the plots. The graphics packages included ggplot2 and factoextra. Biosignal processing and analysis were conducted using MatLab v. R2021b from The Math-Works Inc.