Digital biomarkers for interstitial glucose prediction in healthy individuals using wearables and machine learning

Study design

Ethical approval for the study was obtained in Germany from the corresponding ethics committee of the University of Luebeck (AZ 21-314 and AZ 2022-550). All participants provided written informed consent in accordance with Good Clinical Practice.

The clinical study consisted of two phases: (1) a main study and (2) a follow-up study (as illustrated in Fig. 6a). (1) For the main study, 32 healthy women and men aged 20–40 were recruited at the University of Luebeck. Notably, the sample size extended the number of participants of prior studies with comparable methodologies^22,34. However, formal sample size estimation was not conducted prior to data collection. Due to the lack of available data sets predicting interstitial glucose levels based on non-invasive wearable data without incorporating secondary parameters such as diet or physical activity, the sample size was determined based on practical considerations. In addition to a comprehensive blood panel and the assessment of various liver, kidney, and thyroid parameters, the levels of hemoglobin A1c (HbA1c) were measured to characterize the healthy cohort and specifically confirm that no (pre-)diabetics participated in the study. The HbA1c level was within the normal range at 5.2 ± 0.2%. Throughout the two-week study phase, subjects wore the Abbott Freestyle Libre 2 sensor, on their upper arm, that measured IG levels (mg/dL) continuously every 15 min. Furthermore, participants self-reported their diet, activity, and sleep with the help of an app ([Perfood, Germany]). Within the two weeks, two study visits took place with a duration of 7–8 h each, in which participants wore two additional sensor devices, (1) a 4-channel BiosignalsPlux explorer kit ([PLUX Biosignals, Lisbon]), that collected four sensor modalities (with a sampling rate of 200 Hz): Body temperature ‘BTEMP’, respiratory rate ‘PZT’, electrooculography ‘EOG’ and electrogastrography ‘EGG’, and (2) an Empatica E4 wristband ([Empatica, Boston), that recorded another 4 sensor modalities including skin temperature ‘STEMP’ (at 4 Hz), electrodermal activity ‘EDA’ (at 4 Hz), photoplethysmography (PPG) (producing blood volume pulse ‘BVP’ (64 Hz) and heart rate ‘HR’ (1 Hz)) (sensor placements are exhibited in Fig. 1b). Two standardized, isocaloric, isovolumetric test meals were consumed per visit in a cross-over design: a mixed meal test (MMT) (Fresubin® 2 kcal) and an oral glucose tolerance test (OGTT) (as shown in Fig. 6c and d). Participants who received the OGTT as the first test meal on their first visit received the MMT after a two-hour intervention phase and at least a two-hour washout phase, followed by another two-hour intervention phase. On their second visit, these participants received the MMT first, analogously. The choice of the first test meal (MMT or OGTT) was balanced across subjects. Before and after the different test meals, several blood, urine, and stool samples were taken for analysis of metabolome and other clinical parameters. In total, more than 1550 IG measurements and corresponding 15-min non-overlapping epochs of wearable sensor data containing the eight sensor modalities were collected. In addition, user inputs (i.e., demographic data, including age, biological sex, and BMI), time-domain inputs (seasonal information and corresponding timestamp on consumption), and food diary (standardized meal intervention) were assessed.

A follow-up study (2) was designed to validate the results of the main study under everyday conditions. Five randomly selected participants (P2b, P3b, P6b, P7b, and P10b) of the main study population participated in this follow-up study. Again, the Abbott Freestyle Libre sensor now with a higher sampling rate (every 5 min) continuously recorded their IG levels over two weeks. Throughout the two-week study period, the participants wore the Empatica E4 wristband and ate according to their individual habitual diets.

In total, over 14,400 IG measurements and corresponding 5-min non-overlapping epochs of wearable sensor data limited now only to five sensor modalities, user inputs and time-domain inputs were measured.

Data preparation

High-quality data preprocessing is essential for robust glucose prediction, as it addresses key challenges like sensor noise, temporal misalignment, and physiological variability. Proper techniques such as artifact removal, adaptive resampling, and multimodal signal synchronization significantly improve model inputs while preserving critical glycemic patterns^60,61,62.

For this study (see Fig. 7), the following preprocessing steps were performed: resampling, denoising, and segmentation. To address temporal inconsistencies in the sensor data, including missing values, irregular sampling intervals, and varying recording frequencies across different devices, we developed a systematic data alignment protocol. The approach involved creating a master timeline with fixed, equidistant intervals determined by both the required feature sampling rate and the prediction horizon. This uniform temporal framework served two critical purposes: first, it provided a standardized reference for synchronizing all sensor data streams; second, it enabled precise temporal mapping of each data point relative to the target glucose measurements. The sensor modality data were then interpolated with a common sampling rate and their respective time information to match the predefined time steps. This method effectively normalized the temporal representation of our heterogeneous dataset while preserving the integrity of inter-signal dynamics essential for accurate glucose prediction.

After reviewing the functional subbands of all modalities, it was determined that all sensor data could be resampled to 20 Hz. Resampling reduces training runtime while maintaining performance for practical glucose prediction. Some sensors, like heart rate (1 Hz) or IG (low sampling rate), require upsampling.

Three types of interpolation approaches, such as Linear Interpolation (LiI), Nearest Neighbors’ Interpolation (NNI), and Lagrange Interpolation (LaI), were tested⁶³. Considering the long PH interpolated, the advantage of LiI in the frequency domain⁶⁴, and the requirement of computational resources, we chose to deploy the LiI method for the resampling task, except for IG. For instance, LiI was employed for sensor modalities recording such as EDA and temperature, as these parameters fall within a numerical scale range. Referring to Eq. (1), LiI is a straightforward technique used to calculate values at points that lie between known data points. Mathematical derivation involves finding a linear polynomial that intersects the unique points \(\left({x}_{0},{y}_{0}\right)\text{ and }({x}_{1},{y}_{1}).\)

where \(\varepsilon_{LI} \left( x \right)\) is commonly known as the first-order polynomial interpolation function. Due to its relatively low computational complexity, LiI is frequently used for estimating intermediate data.

In contrast, the nearest neighbor interpolation (NNI) method was tested (refers to Eq. 2) to estimate the value of a new point based on the value of the closest existing point in the dataset. Unlike more complex methods, it does not involve calculating intermediate values or derivatives. Instead, it simply finds the nearest data point and assigns its value to the new point.

Here, \(\left( {xi,yi} \right)\) is a set of discrete datapoint with \(i = \, 0, \, 1, \, 2, \ldots ,n\). \(\varepsilon_{NN} \left( x \right)\) represents the interpolated value at any point \(x\). \(argmin_{i} \left| {x – x_{i} } \right|\) denotes the index k for which the distance \(|x-{x}_{i}|\) is the smallest.

Notably, LiI lacks smoothness at subinterval boundaries due to its non-differentiability at these points. NNI, while simple and fast, creates a piecewise constant function that can result in significant errors between data points with large value differences. To avoid these shortcomings, for glucose-labeled data interpolated using PH as the distance between two consecutive data points, we tested Lagrange Interpolation (LaI) (see Eqs. 3 and 4), in addition to the LiI and NNI methods outlined above. Specifically, LaI generalizes linear interpolation by using higher-order polynomials (the polynomial is set to order five). For a set of n + 1 data points, the Lagrange polynomial is formulated as:

where \(\beta_{i} \left( x \right)\) is the Lagrange basis polynomial defined as:

In addition, body temperature (BTEMP) was calculated accompanied by skin (STEMP) and core temperatures (CTEMP) refers to a proposed and verified criteria (Eq. 5) by R. Lenhardt et al.³⁸ with a predefined \(\alpha\) = 0.67 since our follow-up study was performed in the summer time (daylight saving time in Europe):

Wearable measurements may also encounter outliers due to various factors. Outliers exceeding four standard deviations from the mean are identified using Eq. (6).

where \(x\) is a single sensor measurement, \(\mu\) and \(\sigma\) represent the mean and standard deviation across all values in a sensor modality for each participant, respectively.

Abnormal outliers in the dataset were replaced with means, using a threshold of four deviations from the mean. This choice was based on the empirical rule that 99.7% of data points lie within three standard deviations, with an additional standard deviation buffer to accommodate possible non-normal distributions in sensor data.

To prepare the data for input into ML models, the temporal information from the time series \(T = \left[ {t_{0} ,t_{1} , \ldots ,t_{i – 1} } \right]\) for all sensor modalities \(S = \left[ {s_{0} ,s_{1} , \ldots ,s_{n – 1} } \right]\), given by the devices and used for the experiment along with the measured data, must be defined to fit within the PH’s temporal length, represented as \(p\) in min, which corresponds to the labels \({\Omega } = \left[ {\omega_{0} ,\omega_{1} , \ldots ,\omega_{m – 1} } \right]\). The sensor data from the wristband, which serves as our training data, is partitioned into equidistant time windows \(J = \left[ {j_{0} ,j_{1} , \ldots ,j_{m – 1} } \right]\), with the temporal length of these windows determined by p. For the label \(\omega_{x}\), the time interval \(\left[ {t_{{\omega_{x – p} }} ,t_{{\omega_{x} }} } \right]\) forms the time window \(j_{x}\) for all \(S\). Depending on the length of the PH (\(p\)), a varying quantity of time windows arises. Consequently, when the time series length T of the label data remains constant, the quantity of time windows J is inversely correlative to the growth of \(p.\) The formula for calculating the number of windows is shown in Eq. (7).

A brief PH might lack information, while an extended one could reduce prediction performance. Therefore, a 15-min PH was chosen, with a 5-min PH for comparison purpose in experiments.

Correlation analysis

After data preprocessing, a basic data-driven correlation ablation test was performed to comprehend the data and explore the feasibility of glucose prediction for typical meals using non-invasive sensor records. Without any FE, we adopted the tree-based (decision tree and random forest) and gradient-boosting (GB)-tree-based approaches, which are widely used and effective in various healthcare research studies^{22,65,66,67,68}, with the R² metric based on 32-fold cross-validation to scrutinize the inter- and intra-relationships between interstitial postprandial glucose alterations and different combinations of eight sensors. To mitigate analysis bias between specific models, decision tree and random forest were used for the tree-based category, while LightGBM and XGBoost were deployed for the GB-tree-based category to calculate the average fit score R². As the appropriate research period for the glucose metabolic response is two hours post-meal, which is precisely the length of our data recording, we segmented them in terms of meal type (MMT and OGTT, excluding the washout period). Afterwards, it collaborated with polynomial-interpolation-based glucose meter values (1 value per minute) to reveal the variances and fit patterns between the given sensor modalities or their combinations and the glucose responses for a given meal type consumed in that 2-h interval. This procedure establishes a baseline for subsequent feature-driven continuous glucose prediction studies.

Feature engineering

We developed a population predictive model based on engineered features, which can be used as a general technique to maximize the benefits of limited datasets and incorporate domain expertise into the machine learning (ML) process, with LOPOCV to predict glucose levels in healthy individuals.

A total of 45 types of features across four domains were extracted (see STable 3): demographics, statistical order, time- and frequency/amplitude-domain features, with a PH of 15 min. Z-scores⁶⁹ were calculated (Eq. 8) to address individual variations, enabling effective prediction and calibration using both z-scores and actual measurements.

where \(\Gamma\) indicates the standard z-score calculated from actual measurements \(x\), mean \(\mu,\) and standard deviation \(\sigma.\)

Forty-five features encompassed demographic data such as age, sex, and BMI, indicating biological values. These variables, with dimensions of 5 × 32 for input case (IC) in the main study and 5 × 5 for follow-up study input case 1 (FUIC1) (STable1), were kept constant per participant. Each 15-min interval contained 14 statistical order features, including mean, maximum, minimum, standard deviation, median, zero crossing, the temperature difference between STEMP and BTEMP, 20th, 50th, and 80th percentile of quartile, the interquartile difference between 25 and 75th percentile, and autocorrelations, calculated for each sensor modality (e.g., with dimensions of 14 × 5 × 32 for IC2 in the main study and 14 × 5 × 5 for FUIC1 in the follow-up study).

Furthermore, advanced signal processing metrics were applied to features such as pit and peak calculation, discrete wavelet coefficients, inverse fast Fourier transform coefficients, kurtosis, spectral energy, and spectral entropy, generating 23 frequency and amplitude domain features. These features were applied uniformly to all sensor data, resulting in dimensions of 23 × 5 × 32 for IC2 and 23 × 5 × 5 for FUIC1.

Our emphasis on time-domain features, specifically ‘minutes_from_midnight’ (current day’s minutes from midnight to encode circadian rhythm) and ‘days_from_2021’ (days from 2021’s start to calculate seasonal rhythms), addresses their significant impact on glucose levels. The former not only signifies meal timing but underscores the circadian rhythms’ physiological relevance. Additionally, ‘days_from_2021’ is crucial for accounting for seasonality. This feature accommodates variations in ambient temperatures, influenced by daylight saving time, impacting data collected by wearable sensors. These factors, alongside clustering patterns in daily life, contribute to understanding the practical implications of these features in predicting and calibrating glucose responses effectively.

As we know, diet has a large impact on glucose responses. From a physiological perspective, including meal information in glucose prediction is logical. However, in daily life, thousands of complex composition diets make quantitative and objective analysis arduous. In addition, self-reported food diaries are complex, time-consuming, and error prone. Therefore, we only implemented the food diary feature ‘consumed_food’ to facilitate the correlation between standardized meals and sensor data under predefined experimental circumstances in the main study. After the correlation was established, we eliminated the interference of food factors in our follow-up study to examine whether the sensor data could realistically reflect the underlying logic of glucose responses and whether only using sensor data could effectively predict glucose values in real-time under the unrestricted conditions of daily life.

PPGR prediction models

In the main study, our proposed model trained and tested more than 1550 IG meter measurements and corresponding 15-min PHs from wearables in combination with a food diary (standardized meal intervention: MMT and OGTT) over a period of approximately 14 h (including washout phase) in two days across 32 participants. Subsequently, a dataset of more than 14,400 glucose meter measurements and corresponding 15-min PHs of non-invasive sensor data without food diaries and activity data measured over ten days across randomly selected five participants was exploited in the follow-up study.

We tested two regression models for predicting 15-min PHs of IG measurements based on 45 hand-crafted features implemented with various libraries (sklearn⁷⁰ and scipy⁷¹) in Python with TensorFlow framework⁷² for main and follow-up studies, i.e., RF (with the configuration of hyperparameters: n_estimators = 100, max_depth = 16, n_jobs = − 1) and LightGBM (with the hyperparameters’ configuration: boosting_type = ‘gbdt’, objective = ‘regression’, num_leaves = 5000, learning_rate = 0.1, n_estimators = 100, max_depth = 16, metric = ‘mse’, bagging_fraction = 0.6, feature_fraction = 1.0, reg_lambda = 0.9, importance_type = ‘gain’). A nonFE-based recurrent neural network called LSTM was used for performance comparison. After analyzing the experimental outcomes, it was noted that the RF model performed marginally better in the main study, which involved a limited dataset. In contrast, the LightGBM predictive model was superior in processing large-scale daily data in the follow-up study with a significantly lower prediction error despite even the absence of food diaries. These experiments provided strong evidence for the effectiveness of the predesigned ensemble FS strategy that integrates the advantages of RFE and Boruta (called BoRFE), with LOPOCV and the strength of engineered features.

Typically, FS is adopted to select the effective feature set based on the contribution of each feature in the algorithm, aiming to increase the algorithm interpretability and reduce the model complexity by identifying influential features. Unfortunately, finding those optimal feature sets, especially for high-dimensional data, is a challenging task, even though many ML models are capable of generating FR. Features have a complex relationship regarding their importance to the target variable or classification, so using the features with the lowest score may not lead to optimal results. Consequently, we proposed BoRFE, which involves a two-layer FS process: Boruta effectively eliminates irrelevant features in large feature spaces. At the same time, RFE further refines the selection by iteratively removing one feature at a time collaborating with LOPOCV, which can remove the need for manual trial and error exploration. In the filtering process of BoRFE, only highly relevant hand-crafted features were always involved in the glucose prediction model’s training and optimization to ensure the experimental outcomes’ high interpretability.

Evaluation metrics

RMSE and MAPE were used as quantitative assessment metrics for evaluating glucose prediction models for each fold of LOPOCV. Meanwhile, CEGA was utilized as a qualitative evaluation to visualize glucose prediction error distributions, which is an essential tool for examining the clinical accuracy of blood glucose self-monitoring. Its space is typically separated into five zones, i.e., A, B, C, D, and E, the brief description is as follows:

• A: Clinically correct treatment decisions.

• B: Clinically benign treatment or no treatment needs.

• C: Overcorrection.

• D: Misdiagnosing between healthy, pre-diabetes, and diabetes patients.

• E: Completely opposite or wrong diagnosis.

To enhance the interpretability of the experimental outcomes, we performed FR using two criteria: MDI and Gain (number of estimators of both models are fixed at 100), associated with the RF and LightGBM algorithms, respectively. The MDI sums and averages the actual decrease in node impurity across all trees in RF (weighted by the number of samples split), rather than counting splits. At the same time, the result of the importance calculation based on the Gain criterion includes the total gains of the splits using the feature. All features extracted from the recorded data can generally be divided into four main categories: ‘demographics’, ‘time-domain’, ‘phase marker’, and ‘sensor data’, and then further divided into 13 functional subcategories in terms of their specific functions: ‘biological age’, ‘biological sex’, ‘biological BMI (height and weight)’, ‘skin temperature’, ‘blood volume pulse’, ‘heart rate’, ‘electrodermal activity’, ‘body temperature’, ‘respiratory rate’, ‘electrogastrography’, ‘electrooculography’, ‘meal study phase’, ‘ consumption time’ and ‘seasons’. In our proposed BoRFE-based prediction model, feature importance across categories is calculated and evaluated by averaging the importance of each fold in LOPOCV. The FR results are presented as a percentage of the total significance (all importance scores are summed to 1).