Real-time forecasting of blood glucose (BG) levels has the potential to drastically improve management of Type 1 Diabetes, a widespread chronic disease affecting the metabolic system. Most notably, if hypo or hyperglycemia episodes (i.e. glycemic excursion below or above a safe range) could be accurately predicted, then the patient could be timely warned, thus enabling proactive countermeasures to avoid these dangerous conditions. In this work, a novel personalized algorithm for the real-time forecasting of BG is developed by combining the output of a shallow feed forward neural network with an error imputation module composed by an ensemble of trees. Past glucose readings as well as insulin, meals and work/sleep time information are carefully handled to train and boost the prediction performance of the algorithm. The root mean square error over the 6 subjects achieves a mean value of 18.69 mg/dL and 32.43 mg/dL for 30- and 60-minute prediction horizon respectively.
Type 1 diabetes (T1D) is a metabolic disease characterized by the
destruction of the pancreatic cells responsible for insulin production
and thus resulting in an impaired Blood Glucose (BG)
homeostasis. Diabetes treatment relies on external insulin injection aimed at
maintaining BG levels within a physiological range [
Measures to mitigate/avoid these conditions include the administration of exogenous insulin to reduce hyperglycemic excursion or the consumption of fast-acting carbohydrates for hypoglycemic events. Unfortunately, both insulin injection and carbohydrates consumption affects BG only after a considerable amount of time: 45 minutes for insulin and at least 15 minutes for carbohydrates.
Therefore, accurate BG prediction would be of paramount
importance to enable timely corrective actions and thus ensure successful
BG control. This holds true both for standard therapy, where
corrective actions are manually performed by patients, and in automated
and semi-automated systems such as an artificial pancreas [
Nevertheless, BG prediction is by no mean a trivial task. BG
concentration is the result of complex non-linear dynamical interaction
of multiple physiologic subsystems and is influenced by several
factors, often hard to measure. They include timing and magnitude
of carbohydrates consumption, protein and fat content of the meal,
length, intensity and even type (aerobic vs anaerobic) of physical
exercise, stress, illness or menstrual cycle. Furthermore, the modeling
of BG dynamics is hindered by the large variability in the
physiological metabolic response of different individuals. For these
reasons, non-linear and personalized models can represent one of the
best options to address the task. The introduction of continuous
glucose monitoring (CGM) sensors in T1D care has encouraged the use
of data-driven models based on past BG readings, whereas the
further availability of data offered by infusion pumps and fitness bands
opened the possibility of using exogenous inputs as additional
features for describing the model [
Over the last two decades, several non-linear algorithms have been
tested in this framework, including support vector machine [
Up to the present, none of the these models has stand out from the others in terms of prediction accuracy. This consideration lead our group to focus more on feature manipulation, selection and hyperparameters optimization and to explore the possibility of combining different kinds of non-linear learners. In conclusion, the aim of our work is to synthesize an accurate model of BG dynamics by starting from a simple, feedforward NN and to investigate: the impact of hyperparameters optimization and feature selection; the improvement achievable by combining the NN with an error imputation module (EIM) based on a regression trees ensemble. 2
Our study is based on the real patient data provided by the
OhioT1DM Dataset, described in detail in [
Some of these signals are quasi-continuously measured, like CGM, acceleration or basal insulin. Some others are impulse-like, e.g. self-monitoring BG or meal consumption, which are provided only a few number of times along the day. Some others, like work intensity or sleep quality, are instead defined for time windows of the order of some hours.
Some of the impulsive-like features, e.g. insulin boli, have an impact on glucose dynamics that could last up for several hours. In order to include this information in the framework of feedforward NNs, which have no memory about the dynamics of inputs, we had to rely on new features. We described the insulin-on-board (IOB) by convolution of insulin boli and basal with a 6 hours activity curve, whereas the convolution of consumed carbohydrates with an absorption curve (different for slow- and fast-acting carbohydrates) returned the carbohydrates-on-board (COB). These two variables carry information about the dynamics of slow insulin absorption and carbohydrates slow impact on BG, hence they are suitable for being used with feedforward NNs. In a similar way, we described the physicalexercise-on-board by low-pass filtering with second order transfer function the physical activity intensity.
We also introduced the slope of CGM, computed by using the last 2 hours of readings, since the trend of the glycemic profile resulted to be significantly correlated to future CGM readings in the training set. Other tested features include daytime, time and amount of last carbohydrates intake and several filtered version of the original signals. 3
As shown in Figure 1, the proposed model is composed by two parts. The first one is a shallow NN, which is the main predictor. It is trained to predict future BG values with a certain prediction horizon (PH). The second predictor is based on an ensemble of trees. It is called error imputation module, since it is trained to predict the error that is committed by the shallow NN. Finally, as shown in Figure 1, the prediction of the proposed algorithm is obtained by combining the output of the shallow neural network CGMs(t + P H) with the output of EIM, e^(t + P H), to have an accurate value for the expected glucose concentration.
The feature set for the shallow NN is manually determined on a
population level, i.e. by looking for a unique set of feature which can be
used by all subjects. The first criterion for feature selection consists
in excluding all those features that presents too many missing
values in the training set and hence cannot be considered reliable. For
instance, COB was discarded because, in some subjects, the
information about meal is missing for a large part of the training set (e.g.
subject 567 reported only 31 of the expected 141 meals). A second
criterion consists in excluding the features that are expected to have
a negligible impact on prediction accuracy, e.g. the state of illness
or work. This selection was performed with domain experts and by
means of some preliminary evaluation on the training set. The
selected features still presented many missing values that could hinder
the training procedure. Thus we performed, exclusively on the
training set, a first order interpolation on any gap of samples shorter than
30 minutes. Finally, we investigated how many past CGM readings
should be used as inputs in order to improve model accuracy. To this
purpose we performed a bayesian optimization [
Similarly, we determined the number of hidden layers and number of neurons in each layer with an exploratory optimization, which was also performed via bayesian optimization and using the same training/validation split we employed for feature selection. The best performances on the validation set were found when using a single hidden layer with a reduced number of inner nodes (only 5).
While the architecture of the net was optimized on a population level and it will be the same for every subject, the weights of the net were trained individually for every subject, meaning that the resulting model is tailored on each of them. The shallow net is trained on the whole original training set, with its target being either the 6- or 12-steps ahead prediction of CGM, i.e. prediction horizons of 30 and 60 minutes, respectively. Inputs and targets of the net are normalized by the mean and standard deviation computed on the training set. 3.2
We noticed that the prediction of the main NN is affected by a large
error when abrupt changes occur in its inputs. For instance,
consecutive CGM readings showing a large difference in values are often
associated with poor BG predictions. However, other explanatory
factors for the prediction error could be found in those features which
were not used by the shallow NN. The key idea behind the EIM is
to provide an estimate e^(t + P H) of the true error, e(t + P H),
affecting the prediction of the shallow net. The first step to build this
module is to create a new pool of feature, named Corrective
Feature Set, containing the first order differences, at several time lags,
of: CGM, IOB, COB, sleep/work period, skin temperature and
acceleration data. Then, to take into account the large inter-individual
variability, a feature selection step (based on ReliefF [
The delay existing between the target signal and the predicted one can be computed as the temporal shift that minimizes the square of the mean quadratic error between these two signals: 3.3
The effectiveness of the proposed approach is assessed by comparing the predicted profiles with the ones obtained by a benchmark predictor: a shallow neural network based on CGM data only (CGM-NN). This network resorts the structure of the main predictor (i.e. a single hidden layer with 5 inner nodes) but it exploits past CGM data only (16 samples) as input features. This model will also be personalized for each patient. 4
The accuracy of the predicted profiles is evaluated by using four metrics. The Mean Absolute Error (MAE) is defined as:
M AE = 1 XN(y(t) N t=1 y^(tjt
P H)) where P H is the prediction horizon, y(t) is the current CGM reading, N is the length of the whole signal y and y^(tjt P H) is the the PH-steps ahead prediction using the information available up to instant t. Similarly we can define the Root Mean Square Error (RMSE): v RM SE = tuu N1 t=1
N X(y(t) y^(tjt
P H))2 and Coefficient of Determination (COD):
COD = 100 (1 jj(y(t) y^(tjt jj(y(t)
P H))jj22 ) y(t))jj22 where y(t) is the average value of y. COD counts for the variance explained by the predictive model with respect to the total variance of the signal. Its maximum value is 100%. delay = arg min j2[0;P H] N h 1 N P H
X ((y^(tjt t=1
P H) + j) y(t))2i: 5
Since the shallow net is initialized with random weights, results will be reported in terms of mean and standard deviation of the various metrics on 10 different initialization of the algorithm.
We compared two models: one is the benchmark, shallow NN using exclusively past and present CGM readings (CGM-NN); the other is the shallow net employing the Predictive Feature Set and the error imputation module (NN-EIM). Table 1 and Table 2 reports the results obtained with these two models, on each subject, for 30 and 60 minutes prediction horizons respectively. The last row of the tables averages the mean values of the metrics on every subject.
As aforementioned in Section 3.1, no kind of operation was performed on the testing set in order to impute missing values. Table 3 reports the number of CGM samples available for each subject and the number of those predicted by our CGM-NN and NN-EIM.
NN-EIM achieves better results on each subject for all the evaluation metrics, both for PH = 30 and PH = 60 minutes. On average, the RMSE improves from 19.50 mg/dL with CGM-NN to 18.63 mg/dL with NN-EIM on the 30 minutes PH (p-value=0.031) and from 34.26 mg/dL with CGM-NN to 32.27 mg/dL with NN-EIM (p-value=0.031). The p-values are computed with a Wilcoxon signed rank test and show that the improvement, albeit small ( 5% in both cases), is statistically significant with 1 = 0:95 confidence level.
Figure 3 show the boxplots and scatter plots of the average RMSE values for every subject, for a PH of 30 minutes (Figure 3a) and 60 minutes (Figure 3b). The color of the lines linking the scatter plots indicates the magnitude of the difference in RMSE between the two strategies: green-shaded lines mean an improved accuracy from CGM-NN to NN-EIM; viceversa, the lines are red if model accuracy worsen. Both for 30 and 60 minutes ahead predictions, the use of exogenous inputs and EIM results in a systematic improvement of performances.
The RMSE of the NN-EIM ranges from 16.50 mg/dL up to 21.29 mg/dL with PH=30 for subject 544 and 584, respectively. This difference might be related to the occurrence of large oscillations in CGM data which cannot be clearly explained by any input variable in the dataset. Considering a PH = 30 min, the lowest improvement in terms of RMSE is for subject 552 (16,81 mg/dL vs 16.51 mg/dL, CGM-NN vs NN-EIM respectively). The largest one is for subject 540 (21.66 mg/dL vs 20.42 mg/dL, CGM-NN vs NN-EIM respectively). One of the main reasons linked to this marginal improvement is the large number of missing values in the testing set for some features (e.g. acceleration and skin temperature).
Nevertheless, whenever the feature set is reliable — as for subject 544 in Figure 2 — the proposed algorithm (green line) reduces the prediction delay as well as the error when CGM data show a positive increment related to any external input (i.e. CHO ingested, in this case). Furthermore, the use of lagged first order differences of past CGM data provides a prediction which is more adherent to the target CGM than the one obtained by CGM-NN. This is confirmed by an enhanced COD for subject 552 (90,72% vs 91,38%, CGM-NN vs NN-EIM and by the delay (20 min vs 17,5 min, CGM-NN vs NNEIM).
Same conclusions can be found by considering a PH = 60 min. As before, the lowest improvement in terms of RMSE is for subject 552 (30,45 mg/dL vs 29,56 mg/dL, CGM-NN vs NN-EIM). The largest is for subject 584 (36,85 mg/dL vs 33,84 mg/dL, CGM-NN vs NNEIM). The prediction capabilities of the proposed approach are also confirmed by COD and delay, even when a marginal improvement in terms of RMSE is found. In fact, for subject 552 we found the COD for CGM only NN vs proposed algorithm is 70,12% vs 72,93%. Delay is reduced at 40 minutes, by meaning that the prediction has a useful time anticipation of 20 minutes.
(a) PH=30 minutes (b) PH=60 minutes In this work we presented a new approach for a real-time forecasting of glucose levels based on a shallow neural network and an error imputation module (NN-EIM). Comparing the performance, at PH = 30 and PH = 60, of the novel algorithm vs CGM-NN, we demonstrated that accurate feature manipulation and selection steps can effectively improve the prediction accuracy and reduce the delay affecting the predicted profile. However, the presence of many missing values in some variables reduce the improvement brought by the proposed approach. Finally, a further development of this work is the investigation of patterns within the CGM time series and the inclusion of such physiological priors into the proposed algorithm could results in improved performance.
This work was partially supported by MIUR, (Italian Minister for Education, under the initiatives “Departments of Excellence” (Law 232/2016) and “SIR: Scientific Independence of young Researchers”, project RBSI14JYM2 “Learn4AP: Patient-Specific Models for an Adaptive, Fault-Tolerant Artificial Pancreas”.
A repository containing the code developed for this work is available at: https://github.com/jp993/BGPC_2020