Long-Term Glucose Forecasting Using a Physiological Model and Deconvolution of the Continuous Glucose Monitoring Signal.

(1) Objective: Blood glucose forecasting in type 1 diabetes (T1D) management is a maturing field with numerous algorithms being published and a few of them having reached the commercialisation stage. However, accurate long-term glucose predictions (e.g., >60 min), which are usually needed in applications such as precision insulin dosing (e.g., an artificial pancreas), still remain a challenge. In this paper, we present a novel glucose forecasting algorithm that is well-suited for long-term prediction horizons. The proposed algorithm is currently being used as the core component of a modular safety system for an insulin dose recommender developed within the EU-funded PEPPER (Patient Empowerment through Predictive PERsonalised decision support) project. (2)
Methods: The proposed blood glucose forecasting algorithm is based on a compartmental composite model of glucose-insulin dynamics, which uses a deconvolution technique applied to the continuous glucose monitoring (CGM) signal for state estimation. In addition to commonly employed inputs by glucose forecasting methods (i.e., CGM data, insulin, carbohydrates), the proposed algorithm allows the optional input of meal absorption information to enhance prediction accuracy. Clinical data corresponding to 10 adult subjects with T1D were used for evaluation purposes. In addition, in silico data obtained with a modified version of the UVa-Padova simulator was used to further evaluate the impact of accounting for meal absorption information on prediction accuracy. Finally, a comparison with two well-established glucose forecasting algorithms, the autoregressive exogenous (ARX) model and the latent variable-based statistical (LVX) model, was carried out. (3)
Results: For prediction horizons beyond 60 min, the performance of the proposed physiological model-based (PM) algorithm is superior to that of the LVX and ARX algorithms. When comparing the performance of PM against the secondly ranked method (ARX) on a 120 min prediction horizon, the percentage improvement on prediction accuracy measured with the root mean square error, A-region of error grid analysis (EGA), and hypoglycaemia prediction calculated by the Matthews correlation coefficient, was 18.8 % , 17.9 % , and 80.9 % , respectively. Although showing a trend towards improvement, the addition of meal absorption information did not provide clinically significant improvements. (4)
Conclusion: The proposed glucose forecasting algorithm is potentially well-suited for T1D management applications which require long-term glucose predictions.

1. Introduction

Type 1 diabetes (T1D) is an autoimmune condition characterized by elevated blood glucose levels due to the lack of endogenous insulin production [1]. People with T1D require exogenous insulin delivery to regulate glucose levels. Current therapies for T1D management require measuring capillary glucose levels several times per day and the administration of insulin by means of multiple daily injections (MDI) or continuous subcutaneous insulin infusion (CSII) with pumps. More recently, the improvement on accuracy of subcutaneous continuous glucose monitoring (CGM) has enabled access to virtually continuous glucose concentrations measurements (e.g., every 5 min), glucose trends, and their retrospective analysis. In addition, real-time CGM devices include alerts and alarms for concentrations outside of specified ranges and/or rapid changes in glucose. Clinical data suggest that CGM can improve overall glucose control, as measured by glycated haemoglobin (HbA1c) [2], and can reduce the burden of extreme glucose values (hypo- and hyperglycaemia) [3]. In addition, CGM technology has opened the door to new technologies for managing glucose levels such as sensor-augmented insulin pumps with low-glucose insulin suspension [4] and the artificial pancreas [5]. One additional benefit of the CGM technology is that it facilitates the forecasting of blood glucose levels, and consequently, enables pre-emptive action to avoid undesired events, such as hypoglycaemia and hyperglycaemia.

Glucose forecasting in T1D is a relatively mature field with several algorithms having been proposed, and comprehensive and extensive reviews published, which provides a taxonomy of the different types of existing algorithms [6,7,8]. In particular, four main types of approaches were identified depending on the type of model being used: physiological models [9], data-based models [10,11,12], hybrid models [13], and control relevant models [14,15,16]. Another distinction that can be done among the existing algorithms is based on the type of inputs being considered. A significant proportion of these algorithms use CGM data as the unique source of information to forecast glucose levels [17], while other algorithms use additional exogenous inputs such as meal intake and insulin doses [9], and only a few of them take into account physical exercise [18]. Furthermore, additional information such as meal absorption information could potentially further improve such accuracy [19]. Regarding the prediction horizon, most of the algorithms focus on short-term predictions (i.e., <60 min) [20], with fewer works focusing on longer prediction horizons [9,21,22]. Of note, it has been shown that for short prediction horizons (e.g., 30 min), differences in the forecasting accuracy among different algorithms is marginal [23]. Although short-term prediction horizons have been proven to be useful for applications such as predictive glucose alerts and low-glucose suspension systems [4], they are usually not sufficient for other applications such as closed-loop insulin delivery [5] and insulin dosing decision support [24]. Hence, there is a need for further investigation of this topic.

Finally, there has been a recent increase in the use of machine learning techniques for glucose forecasting [12,22,25]. However, these techniques tend to require a significant amount of data in order to be trained, which might be a limiting factor in clinical practice. A proof of this interest is the recent organisation of the Blood Glucose Level Prediction Challenge as part of the 3rd International Workshop on Knowledge Discovery in Healthcare Data [26], and a recent publication of a literature review [8].

As a result of this research effort, commercial applications have started to appear embedded within sensor-augmented insulin pumps (e.g., Medtronic MiniMed 640G with Smart Guard; Tandem t:slim X2™ with Basal-IQ Technology) which have been proven to reduce nocturnal hypoglycaemia by using predictive glucose alerts and a predictive low-glucose insulin suspension system [27,28].

In this paper, we hypothesise that a glucose prediction algorithm based on a physiological model of glucose–insulin dynamics and a deconvolution technique using CGM data to estimate the model states is a good candidate for performing long-term glucose predictions. In addition, we hypothesise that adding information about meal absorption can enhance prediction accuracy. The presented algorithm is currently being clinically evaluated as part of the safety system of a mobile-based decision support system for T1D management within the framework of the European project PEPPER (Patient Empowerment through Predictive PERsonalised decision support) [29].

2. Methods

The proposed glucose prediction algorithm is based on the composite minimal model of glucose–insulin regulation in T1D [30] that uses deconvolution of the CGM signal to estimate some of the model states. In particular, the states of the gastrointestinal model are estimated using the technique proposed by Herrero et al., which has been proven to be an effective way to estimate the glucose rate of appearance from mixed meals [31]. Then, meal information (i.e., carbohydrate amount and absorption type) and insulin boluses are considered as exogenous inputs. Note that, although more complex models are available in the literature [32], the employed composite minimal model was chosen because of its trade-off between simplicity, which facilitates the parameter identification task, and ability to capture the complex glucose–insulin dynamics. The effectiveness of the employed composite model was evaluated by Gillis et al. for glucose prediction using a Kalman filter technique [33] and by Herrero and associates for detecting faults in insulin pump therapy [30].

2.1. Composite Minimal Model

The employed composite model of glucose regulation in T1D is composed of the minimal model of glucose disappearance proposed by Bergman and colleagues [34], and the insulin and carbohydrate absorption models proposed by Hovorka et al. [35].

2.1.1. Minimal Model of Glucose Disappearance

The minimal model of glucose disappearance [34] is described by the equations: where is the glucose concentration at time t; is the insulin action; is the glucose rate of appearance from ingested meals; is the plasma insulin concentration; is the fractional glucose effectiveness; is the insulin sensitivity, which models the ratio between endogenous glucose production and glucose uptake; is the insulin action rate; is the basal glucose; V is the distribution volume; and W is the subject’s body weight.

2.1.2. Insulin Absorption Model

The plasma insulin concentration is estimated by means of the subcutaneous insulin absorption model proposed by Hovorka et al. [35], which is described by the following equations. where and are the subcutaneous short-acting insulin compartments, denotes the plasma insulin concentration, the input represents the subcutaneous insulin infusion, is the time to maximum insulin absorption, is the distribution volume of insulin and is the decay rate.

2.1.3. Glucose Absorption Model

The glucose rate of appearance () is calculated according to the gastrointestinal absorption model by Hovorka et al. [35], which is represented by the equations: where denotes the glucose appearance in the first compartment, represents the rate of glucose appearance, the model input denotes the carbohydrate intake amount, is the time to maximum glucose rate of appearance and is the carbohydrate bioavailability.

2.2. Accounting for Meal Absorption

Meal composition has a profound effect on blood glucose levels [19]. Therefore, taking this information into account can potentially enhance glucose estimation and forecasting accuracy. To account for this information in a practical way from the user’s perspective, meals were classified as fast, medium and slow absorption. In particular, fast-absorption meals were considered to have more than of the area under the curve (AUC) of the rate of glucose appearance () profile and appeared within the first two hours after the meal ingestion; slow-absorption meals with less than of AUC of the profile appeared within four hours, and medium-absorption meals in the middle. To take meal absorption information into account within the employed glucose absorption model, the time-to-maximum absorption rate is redefined for the duration of each meal as follows: where is the time-to-maximum absorption rate accounting for the meal absorption; and represent the time shift on the time-to-maximum absorption rate due to different meal absorption rates; and is the absorption type of the meal (fast, medium and slow). In particular, and were fixed to 20 min based on experiments carried out using the UVa-Padova simulator [36]. Figure 1 shows the average profiles corresponding to the fast, slow and medium meals, as per definition above, from the UVa-Padova simulator for a 60 g intake of carbohydrates.

Figure 1

Average profiles corresponding to the fast, slow and medium meals from the UVa-Padova simulator for a 60 g intake of carbohydrates.

2.3. Glucose Prediction Algorithm

The proposed glucose prediction algorithm uses a discretised version of the presented composite model (Equations (1)–(7)). For this purpose, a forward Euler’s configuration with 1 min step size is used to simulate the model.

Let: be the system equations representing a discretised version of the described composite model, where k denotes the sampling instant. Let represent the model states; represent the model parameters; and represent the model inputs, where denotes the amount of ingested carbohydrates, denotes the insulin boluses (units), and denotes the meal absorption type (slow, medium, fast).

First, the model states are initialised as , being the current CGM measurement. Then, Equation (9) is evaluated M times, M being the number of minutes elapsed between two consecutive CGM measurements (e.g., 5 min → five times). In parallel, the model states of the gastrointestinal model (Equations (6) and (7)) are estimated in real-time by applying a deconvolution of the CGM signal using the technique proposed by Herrero et al. [31]. In particular, the glucose rate of appearance () in the second compartment is estimated as: where is the derivative of the glucose measurements calculated as the slope of the linear regression of three consecutive glucose values, and is the insulin action calculated with Equation (2). In order to reduce the influence of the measurement disturbance, the derivative is bounded by mg/dL per min. To further reduce the effect of sensor noise on , a moving average filter is applied to the signal: where is the filtered signal and n is the length of the moving window ().

Then, the glucose appearance in the first compartment is estimated by: where is the derivative of approximated by calculating the slope of the linear regression of three consecutive values.

Then, the estimated model states being used for forecasting purposes are calculated by means of a weighted average between the calculated states with the discretised model ( and ) and the states estimated by deconvolution ( and ) as follows: where is a tuning parameter that allows putting more weight on the model estimation or on the estimation by the deconvolution technique.

Similarly, the blood glucose state being used for forecasting purposes is calculated as: where is a tuning parameter that allows putting more weight on the model estimation () or on the CGM measurement (). Note that and are constant. In particular, we used when identifying the model parameters and when testing the model since these were the values leading to better population results.

Finally, the model states are updated as: and the model is valuated over the predefined prediction horizon (PH) to obtain the forecasted states , where is the desired forecasted glucose value. Figure 2 shows a graphical representation of the described glucose forecasting algorithm.

Figure 2

Block diagram corresponding to the proposed glucose forecasting algorithm. The whole diagram is executed every time a glucose value () (continuous glucose monitoring (CGM)) is received. Then, the physiological model represented by the green blocks is evaluated over the prediction horizon (PH) to obtain the forecasted glucose ().

2.4. Model Parameter Identification

The proposed prediction algorithm can be individualised by identifying some of the model parameters by using retrospective data. Since identification of all model parameters is difficult due to identifiability problems, some of the parameters, which are known to have less inter-subject variability, were fixed to mean populations values (i.e., , V, , , , ) [37], while others were set by using a priori known information from the subjects, such as body weight (W) and basal glucose (). Finally, parameters , and were identified by using an optimisation technique that minimises the mean absolute relative difference (MARD) between the predicted glucose () and the corresponding glucose measurements (). Matlab fmincon constrained optimisation routine was employed for this purpose. Constraints for the identified parameters were min per U/mL, min and min. Note that parameters , and were identified for each one of the evaluated prediction horizons. Table 1 shows the employed values for the model parameters indicating which ones are a priori known and which ones are identified.

Table 1

Values of the parameters used in the forecasting algorithm. * indicates parameters that are identified and indicates parameters that are known from a priori information from the subjects. The rest of the parameters are fixed to mean population values obtained from the scientific literature [31,35].

Parameter	S_G	S_I	G_b	V	V_i	W	t_maxI
Value	0.02	*	∗∗	0.9	1.2	∗∗	*
Units	min^−1	min^−1perμU/mL	mg/dL	dL/kg^2	mL/kg	kg	min
Parameter	t_maxG	k_e	p_2	A_g	Q_1	Q_2	PH
Value	*	1.5	0.02	0.85	0.5[0.7]	0.5[0.7]	30
Units	min	min^−1	min^−1	–	–	–	min

2.5. Baseline Algorithms

In order to compare the performance of the proposed algorithm, validated and commonly employed glucose forecasting algorithms from the literature were chosen. In particular, a third-order auto regression with exogenous input (ARX) method was chosen because it has shown good performance when compared to other algorithms and it is easy to implement [23]. Then, a latent variable with exogenous input (LVX) method was selected because it has showed superiority when compared against other existing techniques and its source code is publicly available [38]. Of note, more recent algorithms, such as [9,21], could have also been chosen, but their complexity of implementation and difficulty in replicating the reported results was the main reason for not doing so.

2.6. Clinical Data Testing

To test the proposed algorithm, a two-week clinical dataset obtained from 10 adult subjects with T1D undergoing a clinical trial carried out at Imperial College Healthcare NHS Trust, London, UK evaluating the benefits of an advanced insulin bolus calculator was employed [39]. In particular, one week worth of data per participant was used for parameter identification purposes (mean ± std: data points) (physiological model-based (PM) algorithm), while the other week was employed for testing purposes (mean ± std: data points). Since no reliable information about meal composition was available, the assumption that breakfast and snacks were fast absorption meals, and lunch and dinner were medium absorption meals, was made. Note that not having such information might limit the benefits of accounting for meal absorption information. In an ideal scenario, CGM measurements are received every 5 min. However, it is common to have gaps in the data. In order to minimise the impact of such gaps in the proposed method, a modified Akima cubic Hermite interpolation method was employed for data imputation purposes.

2.7. In Silico Testing

In order to evaluate the potential benefit of using meal absorption within the proposed algorithm, a modified version of the UVa-Padova T1DM simulator (v3.2) [36] was employed. Intra-day variability was emulated by modifying some of the parameters of the model described in [32]. In particular, meal variability was emulated by introducing meal-size variability (), meal-time variability () and uncertainty in the carbohydrate estimation (uniform distribution between and ) [40]. Meal absorption rate () and carbohydrate bioavailability (f) were considered to vary by and respectively. To account for variability in meal composition, the 33 available meals in the simulator were considered. Note that each cohort (adults, adolescent and children) has 11 different meals. In addition, 16 new meals obtained by using the technique for estimating the rate of glucose appearance proposed by Herrero et al. in [31] were included. Details about the mixed meal library are provided in Appendix A. By using the absorption classification criteria introduced in Section 2.2, of the 49 considered meals, 31 were classified as fast absorption, 15 as medium absorption and three as slow absorption. Finally, intra-subject variability in insulin absorption model parameter (, , ) was assumed as [41,42]. The 10 available adult subjects were used for this purpose. The open-loop insulin therapy provided by the simulator was employed to generate the datasets. A two-week scenario with a daily pattern of carbohydrate dose intake of 7 am ( min) (70 g), 13 pm ( min) (100 g) and 7 pm ( min) (80 g) was chosen. The first week of data was used for the parameter estimation (2016 data points), while the second week works for testing purposes (2016 data points).

2.8. Evaluation Metrics

To measure the forecasting accuracy of the evaluated algorithms, the root mean square error (RMSE), and error grid analysis (EGA) [43], were used. RMSE is calculated as: where is the forecasted glucose value as indicated in Section 2.3, is the corresponding CGM measurement, and N is the total number of glucose data pairs.

EGA is employed to represent the clinical significance of the error between the forecasted glucose value and the actual measurement. In particular, the A region of EGA corresponds to accurate glucose predictions, meaning that forecasted glucose values deviates less than from the actual CGM measurements (i.e., ), or when both the forecasted and the actual measurements indicate hypoglycaemia (i.e., mg/dL and mg/dL). The B region corresponds to forecasted values with a deviation from the CGM measurements larger than (i.e., ) but would not lead to inappropriate treatment (see regions C,D,E). The C region represents forecasted glucose values are more than 100 mg/dL below the actual measurement (i.e., mg/dL), or forecasted glucose is indicating hypoglycaemia while the CGM measurements is between 130 and 180 mg/dL (i.e., mg/dL and 130 mg/dL mg/dL). Finally, the D region represents the zone with forecasted glucose values in the target range while the actual measurement is out of the target range (i.e., 70 mg/dL mg/dL and mg/dL mg/dL). The E region corresponds to potentially clinically dangerous glucose predictions, meaning forecasted glucose indicating hypoglycaemia while CGM measurements indicating hyperglycaemia (i.e., mg/dL and mg/dL), or when forecasted glucose indicate hyperglycaemia while CGM measurements indicate hypoglycaemia (i.e., mg/dL and mg/dL).

In addition, the efficiency of the algorithm at predicting hypoglycaemia was evaluated by the Matthews correlation coefficient (MCC), defined as: where denotes the number of true positives, i.e., predictions of hypoglycemic events that are confirmed to be actual episodes of hypoglycemia, with a hypoglycemic event being defined as three consecutive glucose values below 70 mg/dL; denotes the number of false negatives, i.e., missed predictions of hypoglycemic events; denotes the number of true negatives, i.e., correct prediction of glucose values above 70 mg/dL; and denotes the number of false positives, i.e., false predictions of hypoglycemic events.

Assessment of statistical significance for between-method differences was performed, with , and confidence levels, using a paired t-test as implemented in Matlab.

Finally, prediction horizons of 30, 60, 90, and 120 min were employed to compare the studied algorithms.

3. Results

3.1. Clinical Data Results

Baseline characteristics of the 10 selected participants (one male and nine female) had a median (interquartile range (IQR)) age of 29.5.5 (25.0–42.7) years, duration of diabetes 16.0 (9.2–24.7) years, BMI 26.0 (24.4–31.0) Kg/m, and HbA1c 53.0 (51.0–59.7) mmol/mol.

The distribution of the identified model parameters for the selected 10-adult cohort is min per U/mL, min and min.

Table 2 shows the results corresponding to the 10 adult cohort for baseline algorithms (LVX and ARX) and the proposed method, without accounting for meal absorption information (PM) and accounting for it (). In order to show the consistency of the results at the individual level, the RMSE results for each of the evaluated subjects are presented in Appendix B (Table A4). The normality of the reported data distributions was evaluated by means of the Shapiro–Wilk test and, in all cases, the null hypothesis for an was not rejected.

Table 2

Results corresponding to the 10-adult population expressed as for the studies algorithms and different PH in minutes. Assessment of statistical significance between adjacent rows is indicated with for , for , and for .

Config-uration	PH =30
	RMSE (mg/dL)	EGA-regions (%)					MCC
		A	B	C	D	E
LVX	20.74	85.56	13.25	0.08^ +	1.08^ +	0.03	0.71^ *
	±2.51	±4.55	±4.02	±0.06	±0.71	±0.02	±0.12
ARX	20.74^ *	85.34	13.22^ +	0.06^ ‡	1.36	0.02^ ‡	0.67
	±2.65	±4.49	±3.80	±0.05	±0.93	±0.02	±0.12
PM	19.03^ *	86.21 *	12.65 *	0.03^ ‡	1.10^ +	0.01	0.72^ *
	±2.28	±5.07	±3.73	±0.03	±0.94	±0.01	±0.12
PM_MA	17.67	87.94	11.11	0.02	0.92	0.01	0.74
	±2.12	±4.61	±3.39	±0.02	0.77	±0.01	±0.11
Config-uration	PH =60
	RMSE (mg/dL)	EGA-regions (%)					MCC
		A	B	C	D	E
LVX	34.93	67.62	29.26	0.45	2.47^ ‡	0.20^ *	0.42^ ‡
	±4.85	±6.30	±5.08	±0.27	±1.53	±0.13	±0.14
ARX	34.52^ ‡	67.11	29.39	0.34^ ‡	3.00	0.16^ ‡	0.39
	±4.82	±6.35	±4.76	±0.22	±2.15	±0.12	±0.15
PM	32.69 *	68.10^ *	28.70^ *	0.20^ +	2.90^ *	0.10^ +	0.43^ +
	±4.17	±8.14	±5.64	±0.15	±2.58	±0.08	±0.12
PM_MA	30.36	70.80	26.32	0.12	2.70	0.06	0.44
	±3.88	±8.05	±5.59	±0.10	±2.48	±0.06	±0.12
Config-uration	PH =90
	RMSE (mg/dL)	EGA-regions (%)					MCC
		A	B	C	D	E
LVX	44.56	56.36	38.57	1.04	3.49	0.54	0.30^ ‡
	±7.76	±6.68	±4.66	±0.65	±2.06	±0.46	±0.14
ARX	43.65^ ‡	56.45	38.25	0.87	4.00	0.43^ ‡	0.27^ +
	±7.24	±6.74	±4.54	±0.65	±2.77	±0.33	±0.14
PM	41.07^ *	58.17^ *	37.09^ *	0.53^ *	3.91^ ‡	0.30^ +	0.38^ *
	±5.65	±7.63	±4.96	±0.25	±2.94	±0.21	±0.10
PM_MA	38.14	61.05	34.73	0.33	3.70	0.19	0.42
	±5.25	±7.69	±4.94	±0.18	±3.91	±0.13	±0.11
Config-uration	PH =120
	RMSE (mg/dL)	EGA-regions (%)					MCC
		A	B	C	D	E
LVX	51.30	48.99	44.22	1.53	4.23	1.03	0.22^ ‡
	±10.26	±6.83	±4.31	±1.05	±2.35	±0.97	±0.12
ARX	49.85^ *	49.34^ *	43.89^ *	1.31^ *	4.70^ ‡	0.76^ +	0.21^ *
	±9.33	±6.72	±4.32	±1.07	±2.95	±0.60	±0.12
PM	43.57^ *	55.25^ *	39.45^ *	0.69^ *	4.23^ ‡	0.38^ +	0.35^ *
	±6.53	±7.46	±4.69	±0.35	±3.04	±0.27	±0.09
PM_MA	40.46	58.19	37.08	0.45	4.04	0.24	0.38
	±6.06	±7.58	±4.78	±0.25	±3.01	±0.17	±0.10

Table A4

Individual RMSE results (in mg/dL) corresponding to the 10-adult clinical cohort for the studied algorithms and different PH in minutes.

Config-uration	PH =30
	# 1	# 2	# 3	# 4	# 5	# 6	# 7	# 8	# 9	# 10	Mean	STD
LVX	19.69	16.14	22.91	21.80	18.34	22.50	20.19	23.95	23.06	18.81	20.74	2.51
ARX	20.11	16.02	22.96	21.66	18.08	22.75	20.10	24.39	22.90	18.48	20.74	2.65
PM	17.80	14.92	20.88	20.02	16.67	21.06	19.95	20.40	21.72	16.90	19.03	2.28
PM_MA	16.52	13.85	19.39	18.59	15.48	19.56	18.52	18.94	20.17	15.69	17.67	2.12
Config-uration	PH =60
	# 1	# 2	# 3	# 4	# 5	# 6	# 7	# 8	# 9	# 10	Mean	STD
LVX	35.70	27.61	35.83	36.73	29.08	39.58	33.11	38.77	42.57	30.35	34.93	4.85
ARX	36.46	27.11	35.34	35.96	28.55	39.42	32.71	39.70	40.31	29.61	34.52	4.82
PM	32.15	25.87	34.99	32.28	28.52	37.49	33.93	34.79	38.77	28.07	32.69	4.17
PM_MA	29.86	24.03	32.49	29.98	26.48	34.81	31.51	32.30	36.01	26.07	30.36	3.88
Config-uration	PH =90
	# 1	# 2	# 3	# 4	# 5	# 6	# 7	# 8	# 9	# 10	Mean	STD
LVX	48.43	35.51	40.26	45.55	34.00	52.17	43.48	51.07	57.44	37.69	44.56	7.76
ARX	48.83	34.77	41.20	44.43	33.15	50.57	42.55	51.84	52.70	36.49	43.65	7.24
PM	42.24	34.24	38.85	39.88	34.23	47.22	43.22	45.29	50.54	35.00	41.07	5.65
PM_MA	39.22	31.80	36.08	37.03	31.79	43.85	40.14	42.06	46.94	32.50	38.14	5.25
Config-uration	PH =120
	# 1	# 2	# 3	# 4	# 5	# 6	# 7	# 8	# 9	# 10	Mean	STD
LVX	58.161	41.06	42.03	49.92	37.62	60.52	52.11	60.86	67.83	42.91	51.30	10.26
ARX	58.07	40.27	43.99	48.88	35.99	57.53	50.61	61.14	61.18	40.87	49.85	9.33
PM	45.94	36.53	38.90	42.32	35.22	50.20	46.14	49.04	54.22	37.18	43.57	6.53
PM_MA	42.67	33.92	36.12	39.30	32.71	46.62	42.85	45.55	50.35	34.53	40.46	6.06

Figure 3 depicts the percentage of improvement, for each of the prediction horizons, of three of the evaluated metrics (RMSE, A-region of EGA, and MCC) when comparing the proposed method versus the ARX model, which ranks second in terms of prediction accuracy.

Figure 3

Average percentage of improvement against prediction horizon for three of the evaluated metrics (RMSE, A-region of EGA, and MCC) when comparing the proposed method versus the model on the 10-adult real cohort.

Figure 4 shows a 24 h period close up for a selected individual showing the prediction results for the three evaluated forecasting methods (LVX, ARX, and PM) with a prediction horizon of 120 min.

Figure 4

Example of 24 h period close up for a representative real individual showing the prediction results for the three evaluated forecasting methods with a prediction horizon of 120 min. The continuous glucose measurements is represented by the dashed black line, the prediction by the proposed PM method is displayed in solid-red line, results for the LVX and ARX methods are showed in dotted green line and dash-dotted blue line respectively. Vertical pink bars indicate carbohydrate intakes (grams) and vertical light blue bars indicate insulin boluses (units).

3.2. In Silico Results

The distributions of the identified model parameters for the selected in silico cohort of 10 adults with T1D, expressed as , are: min per U/mL, min and min.

Table 3 presents the results corresponding to the 10-adult cohort for each one of the selected baseline algorithms (LVX and ARX) and the proposed method, without accounting for meal absorption information (PM), and accounting for it (PM).

Table 3

Results corresponding to the 10-adult virtual population expressed as for the studies algorithms and different PH in minutes. Assessment of statistical significance between adjacent rows is indicated with for , for and for .

Config-uration	PH =30
	RMSE (mg/dL)	EGA-regions (%)					MCC
		A	B	C	D	E
LVX	12.85	94.93^ *	4.62^ *	0	0.45^ *	0	0.94
	±1.47	±0.90	±0.81	±0	±0.31	±0	±0.03
ARX	12.84^ *	93.80^ +	5.39^ *	0.01	0.80^ *	0	0.95^2
	±1.31	±1.08	±0.96	±0.02	±0.34	±0	±0.02
PM	10.90^ *	94.82^ *	4.61^ *	0	0.57^ *	0	0.98^ *
	±1.23	±1.19	±0.79	±0	±0.15	±0	±0.06
PM_MA	10.06	95.70	3.92	0	0.38	0	0.99
	±1.14	±1.16	±0.82	±0	±0.21	±0	±0.06
Config-uration	PH =60
	RMSE (mg/dL)	EGA-regions (%)					MCC
		A	B	C	D	E
LVX	28.00^ ‡	74.07	24.26^ ‡	0.25	1.41^ *	0.01	0.79^ +
	±3.62	±4.87	±5.19	±0.37	±0.70	±0.02	±0.08
ARX	26.38^ +	75.91	20.58	0.07	3.43	0.01	0.69^ *
	±3.32	±3.84	±3.54	±0.11	±1.02	±0.02	±0.12
PM	24.44^ *	76.46^ *	20.29^ *	0	3.25^ *	0	0.83^ *
	±2.61	±4.24	±3.37	±0	±0.61	±0	±0.12
PM_MA	22.56	79.13	18.02	0	2.85	0	0.87
	±2.41	±4.00	±3.17	±0	±0.49	±0	±0.12
Config-uration	PH =90
	RMSE (mg/dL)	EGA-regions (%)					MCC
		A	B	C	D	E
LVX	41.04^ ‡	55.42^ ‡	40.92^ +	1.34^ +	2.20^ *	0.12	0.53^ *
	±6.21	±11.35	±11.93	±0.75	±1.05	±0.10	±0.13
ARX	36.10^ +	64.24	29.85	0.34^ ‡	5.35	0.22	0.47^ *
	±4.92	±5.15	±4.83	±0.33	±1.50	±0.30	±0.15
PM	33.50^ *	64.86^ *	29.97^ *	0.10	4.99 *	0.08^ ‡	0.55^ *
	±3.51	±5.49	±5.38	±0.16	±0.98	±0.08	±0.16
PM_MA	30.93	67.75	27.39	0.05	4.79	0.02	0.58
	±3.51	±5.36	±5.17	±0.10	±0.87	±0.05	±0.16
Config-uration	PH =120
	RMSE (mg/dL)	EGA-regions (%)					MCC
		A	B	C	D	E
LVX	48.50	46.61	48.63^ ‡	2.00^ *	2.41^ *	0.35	0.38^ *
	±8.97	±14.29	±14.69	±0.82	±1.37	±0.33	±0.16
ARX	43.03^ *	55.67^ +	36.13^ ‡	0.75^ +	6.80^ +	0.65^ ‡	0.29^ *
	±6.06	±6.07	±5.34	±0.46	±1.97	±0.62	±0.15
PM	37.63^ *	60.72^ *	33.09^ *	0.17^ ‡	5.87	0.15^ ‡	0.45^ *
	±4.70	±6.62	54.01^ *	54.01^ *	54.01^ *	54.01^ *	±0.15
PM_MA	34.74	63.64	30.52	0.07	5.69	0.07	0.51
	±4.34	±5.52	±5.42	±0.16	±1.02	±0.10	±0.17

4. Discussion

Results obtained with both clinical data (Table 2) and in silico data (Table 3) are consistent at showing that, when compared to commonly employed glucose forecasting algorithms (LVX and ARX), the proposed PM algorithm is superior on all evaluated prediction horizons. As observed in Figure 3, such superiority is more evident at longer prediction horizons, and in particular, when looking at the MCC. In the context of predictive glucose alarms, this improvement would be translated into a system with less false alarms and missed alarms, while in an automatic insulin delivering system, it could be translated into better glucose control.

This improvement in forecasting accuracy at longer horizons might be attributable to the use of a physiological model, which is able to better account for the long-term glucose–insulin dynamics when compared to purely data-driven based approaches.

As observed in Table 2 and Table 3, the inclusion of meal absorption information into the PM method did not yield significant improvements in the forecasting accuracy. Although additional meals were included into the UVA-Padova simulator to increase its inter-meal variability, such variability might still not be enough to show the benefits of the proposed approach. Note that the presented mixed-meal model library in Appendix A is per se a scientific outcome that can be used in other applications [44].

It is worth noting that, when comparing against the baseline algorithms, the proposed method shows better performance on the real cohort than on the in silico. Our hypothesis to explain this fact is that our method is better at capturing the complexity of the real world glucose–insulin dynamics as a result of using a deconvolution technique for state estimation.

Although the proposed approach has been proven to be superior to some of existing forecasting algorithms, there are still some scenarios were there is a significant mismatch between the predicted glucose and the actual glucose, e.g., Figure 4 around 15:00 h. This might be explained by the presence of an unaccounted perturbation such as physical exercise. Note that for this particular scenario, a closed-loop control algorithm, such as Model Predictive Control, could lead to sub-optimal results, but this is a limitation inherent to this type of feedback controller.

Finally, unlike recently proposed machine learning methods, which require a significant amount of data to be trained, the proposed methods requires a relatively small dataset (e.g., one week), which makes it practical for clinical applications.

5. Conclusions and Future Work

The proposed glucose forecasting algorithm based on a physiological model and a state estimation deconvolution technique is a potential solution for diabetes technology solutions that require long-term glucose predictions (e.g., 120 min), such as closed-loop insulin delivery and insulin dosing decision support. Future work to further improve the accuracy of the proposed algorithm includes accounting for circadian variations on insulin sensitivity [45] and the short- and long-term effect of exercise on insulin requirements [46]. Finally, online parameter estimation is another potential avenue for investigation.

Table A1

Mixed meals information and bibliographic references.

Meal	Ingredients	Weight	CHO	CHO; Prot.; Fat	Reference
#		(Kg)	(g)	(% energy)
1	Scrambled eggs, Canadian bacon, gelatin (Jell-O)	77	77	45;15;40	[37]
2	White bread, low-fat cheese, sucrose, oil, butter	82.3	111	55;15;30	[47]
3	Fat milk, white rice, low-fat cheese, fructose, pear, bran-cookies, oil	82.3	112.3	55;15;30	[47]
4	Pasta, oil (low fat)	57	75	80;15.4;4.6	[48]
5	Pasta, oil (medium fat)	57	75	56;10.8;33.2	[48]
6	Pasta, oil (high fat)	57	75	37.4;7.2;55.4	[48]
7	Rice, pudding, sugar and cinnamon	65^ 🟉	50.5	74.6;14.2;11.2	[49]
8	Toast, honey, ham, curd cheese, orange juice	65^ 🟉	50.2	26.2;16.5;56.7	[49]
9	Pear barley	59.8	50	79;15;5	[50]
10	Instant mashed potato	59.8	50	78;5.5;4.5	[50]
11	2 slices of bread, 1 and 12 eggs, 1 tea spoon of margarine and orange juice	65	50	49;22.3;28.7	[51]
12	Cereal, coconut, chocolate, fruit and whipping cream	76^ 🟉	93	18;16;66	[52]
13	Oats, coconut, almonds, raisins, honey, sunflower oil, banana, double cream and milk	61.9	121.1	48.6;6.9;48	[53]
14	Same as meal 13	61.9	70.3	28.2;6.6;65.2	[53]
15	Same as meal 13	61.9	50	20;6.1;73.9	[53]
16	Oat loop cereal, milk, white bread, margarine, strawberry jam, orange juice	67^ 🟉	68.8	57;19;24	[54]

🟉 Estimated from Body Mass Index (BMI); CHO denotes Carbohydrates; Prot. denotes Protein.

Table A2

Gastrointestinal model parameters corresponding to the 16 selected mixed meals. Coefficient of variation (%) provided by the Matlab lsqnonlin optimization routine is reported in brackets.

Meal #	k_min	k_max	k_abs	b	d
1	0.0123(1.0)	0.0575(2.7)	0.0388(4.5)	0.6947(1.0)	0.0145(5.3)
2	0.0128(1.0)	0.0176(1.1)	1.4807(4.8)	0.9905(0.2)	0.4388(3.2)
3	0.0110(1.5)	0.0207(4.8)	0.0779(13.2)	0.9593(0.7)	0.1633(4.4)
4	0.0096(1.3)	0.0155(1.5)	0.0719(6.0)	0.7834(0.9)	0.3303(3.2)
5	0.0025(44.5)	0.0128(0.9)	0.1276(8.7)	0.7729(0.5)	0.6081(4.7)
6	0.0675(14.3)	0.0256(13.5)	0.0332(25.4)	0.3917(4.6)	0.7238(3.5)
7	0.2000(28.5)	0.0396(2.2)	0.0104(1.8)	0.4048(7.1)	0.4955(2.4)
8	0.0129(3.9)	0.0224(22.2)	0.0273(45.7)	0.7375(8.5)	0.1964(27.8)
9	0.0098(1.7)	0.0335(5.9)	0.0509(11.8)	0.8031(0.9)	0.1979(3.2)
10	0.0326(2.8)	0.1422(6.5)	0.0355(3.8)	0.9003(1.1)	0.0547(7.5)
11	0.0199(1.8)	0.0896(7.6)	0.0374(5.5)	0.9224(0.6)	0.0939(5.7)
12	0.0115(0.2)	0.0200(0.6)	1.6872(1.5)	0.8578(0.4)	0.3289(1.5)
13	0.0138(0.7)	0.0201(0.9)	0.0946(3.6)	0.8443(0.5)	0.3668(2.3)
14	0.0102(0.4)	0.0211(0.6)	1.6260(1.6)	0.9939(0.04)	0.3283(0.6)
15	0.0110(0.8)	0.0196(0.5)	0.1563(4.8)	0.9982(0.02)	0.4613(1.0)
16	0.0104(2.3)	0.0215(1.0)	2.9750(3.9)	0.9028(0.4)	0.5045(2.3)

Table A3

Metrics to evaluate the model fitting to reference profiles.

Meal	△Ra_peak	△T_peak	△AUC	RMSE	R^2
	mg · min^−1· kg^−1	min	%	mg · min^−1· kg^−1	-
1	0.16	1	2	0.1521	0.995
2	0.27	7	1.6	0.31623	0.978
3	0.20	12	1.6	0.3458	0.952
4	0.14	10	1.2	0.17675	0.979
5	0.04	12	3.4	0.2535	0.964
6	0.20	10	12.5	0.37702	0.818
7	0.22	6	3.6	0.230776	0.970
8	0.25	16	4.3	0.17117	0.972
9	0.22	5	5.1	0.26335	0.918
10	0.17	5	1.1	0.2928	0.991
11	0.01	3	4.0	0.25221	0.987
12	0.21	2	2.0	0.28166	0.984
13	0.21	9	0.5	0.27094	0.993
14	0.19	8	1.3	0.11685	0.995
15	0.12	6	2.5	0.09367	0.995
16	0.07	1	0.02	0.3495	0.969