Prediction of Blood Glucose Levels for People with Type 1 Diabetes using Latent-Variable-based Model

Abstract

Regulation of blood glucose concentrations (BGCs) is a tough burden for people living with type 1 diabetes mellitus (T1DM). People with T1DM must administer exogenous insulin to maintain their BGC within an euglycemic range. Hyperglycemia (high BGC) and hypoglycemia (low BGC) can occur because of poor BGC management. Using recursively identified models to predict the future BGC values opens novel possibilities for improving the BGC regulation performance by adjusting the dose of insulin infusion, taking rescue carbohydrates, or both. The BGC prediction model can also benefit the development of artificial pancreas systems. In this paper, a latent variable (LV) based multivariable statistical modeling approach is applied to model BGC dynamics and forecast future BGCs. The LV-based model is a powerful linear method to build an empirical model by using the collected data. The model is evaluated with the Ohio T1DM dataset that contains BGCs from continuous glucose monitoring (CGM) sensors, basal and bolus insulin information from insulin pump, and additional information from wristbands or reported by the subject. The results indicate that the LV-based model can predict future BGC values with high accuracy for prediction horizons of 30 and 60 minutes.

-

People with type 1 diabetes mellitus (T1DM), an immune disorder where the pancreas does not produce insulin, must administer exogenous insulin either by injections several times every day or infusion with an insulin pump to maintain their blood glucose concentrations (BGCs) within a safe range (70-180 mg/dL) [ 23 ]. Without effective regulation, people with T1DM may suffer from several long-term complications caused by hyperglycemia (high BGC), such as kidney failure, blindness, and the deterioration of cardiovascular health. They can also have low BGC (hypoglycemia), which may cause dizziness, diabetes coma, or even death because of the lack of energy for the brain [ 4 ].

Artificial pancreas (AP) systems are proposed to provide more reliable BGC management by automatically calculating the dose of insulin to infuse by using a closed-loop controller incorporating BGCs measured by continuous glucose monitoring (CGM) sensors, historical infused insulin data, and other available information, such as meal carbohydrates (CHOs) and exercise information [ 15, 9, 16, 18, 2, 10 ]. The closed-loop controller is the critical component in developing an AP system. Model-based controllers have become the preferred option in recent in silico and clinical studies where the future BGCs predicted by a model, historical CGM measurements, administered insulin, and constraint conditions are taken into account when computing the future insulin doses to be infused. The data-driven modeling techniques have been evaluated in many studies because of their computational tractability and identification efficiency.

The empirical modeling technologies for T1DM include linear and nonlinear methods. For nonlinear models, artificial neural networks (ANN) [ 1 ], convolutional neural networks (CNN) [ 26 ], recurrent neural networks (RNN) [ 13, 3 ], and other machine learning and deep learning techniques [ 17, 11 ] have been used to model the glucose dynamics. However, the nonlinear models can only be trained with a large data set, and it can be difficult to develop personalized models or to update the model parameters or structure. For the linear modeling methods, autoregressive model with exogenous inputs (ARX), and autoregressive moving average model with exogenous inputs (ARMAX) [ 7, 19, 25 ] are used to build personalized glucose prediction models with recursively updated model parameters where the future BGC values are modeled as a linear combination of historical measurements from sensors, administered insulin, and other information such as meal CHOs and exercise. Statistical methods based on latent variables (LVs) are demonstrated to have powerful abilities for various data analysis tasks, modeling, and process monitoring [ 20, 21 ], and has been proven as a good alternative linear model for type 1 diabetes [ 24 ].

In this paper, a novel multivariate statistical method proposed by Nelson and MacGregor [ 14, 8 ] where the score vector is estimated using missing data imputation technique is applied to model the glucose dynamics based on LVs derived from principal component analysis (PCA). There are three key steps in developing a BGC prediction model based on the LVs technique. First, a PCA is performed on the gathered data to decompose the data into a linear combination of scores and loadings. Then, the unobserved variables are estimated as conditional mean values computed from the gathered data and new measurements. Finally, the score of new observed data is estimated using incomplete observations and the future BGC values are predicted. The rest of this paper is organized as follows: Section 2 summarizes the pre-processing of the data set. The LV-based method is described and a glucose prediction model is developed in Section 3. The Ohio T1DM data set [ 12 ] is used to assess the performance of the model and the results are given and discussed in Section 4. Finally, some concluding remarks are provided in Section 5. 2 2.1

Data Pre-processing

Data

The Ohio T1DM dataset provided by the Blood Glucose Level Prediction (BGLP) challenge records eight weeks of data collected from six T1DM patients with subject ID: 540, 544, 552, 567, 584, and 596. The data set contains BGC values measured by CGM sensors with a sampling time of 5 minutes, basal and bolus insulin information from the insulin pump, information collected from wristbands, and events (i.e., meal, exercise, work, illness, etc.) recorded by the subjects themselves (see [ 12 ] for details). However, there is no evidence on the accuracy of the subject-reported information, which may degrade the reliability of glucose prediction model and the wristband did not work continuously for a long period of time due to its limited battery power, which causes lots of missed data. Thus, only data collected from CGM sensors and insulin pumps is used in this study.

The insulin infused with an insulin pump includes basal insulin, given continuously since the defined start times, and bolus insulin, which is usually given to mitigate the effects of the rapid raise of BGC around meal times. Basal insulin is recorded as “basal” which is the basal insulin infusion rate and “temp basal” which defines the basal insulin infusion rate in specific time intervals to achieve better regulation of BGCs. There are three types of bolus insulin defined in the dataset: normal, normal dual and square dual. For “normal” bolus insulin, a certain dose of insulin is infused immediately to the patient, while for “normal dual” bolus insulin, a certain percentage of bolus is given to the subject up front and the remaining insulin amount is given gradually over a longer time duration. In this study, half of the dose is supposed to be given to the patient immediately and the remaining half is infused over the following 30 minutes. For “square dual” bolus insulin, the bolus insulin is administered continuously in the given time interval. The bolus and basal insulin values are converted to U=min and aligned to match the CGM sampling time.

The insulin on board (IOB), which represents the insulin that remains active within the body, is calculated from a physiological model [ 22 ] and is found useful in some previous works [ 1, 6 ]. The physiological model is represented as dC1 (t)

dt dC2 (t) dt = u (t)

KdiaC1 (t) = Kdia (C1 (t)

C2 (t)) IOB (t) = C1 (t) + C2 (t) (1) (2) (3) where C1 and C2 define two insulin compartments, u (t) is the insulin infusion, and Kdia is a constant related to the duration of insulin action, set as 0.0182.

The dataset contains missing CGM measurements for in both the training and testing sets that cause discontinuous CGM curves. The gaps that are not greater than 3 samples in the training dataset are interpolated using first-order linear model. For the cases where more than 3 samples are missed, a single variable statistical model similar to the one used to modeling glucose dynamics is used. In the testing data set, with the intent of on-line implementation, only the historical CGM measurements are used to fill the missing gaps until new the CGM observations become available. Fig. 1 shows an example of the filled gap in the training dataset (Subject ID 540). 2.2

Batch Generation

The pre-processed CGM readings (y) and calculated IOB (IOB) are arranged as time series and a sliding window of 2 hours length is chosen to generate batch segments from the training data as 2 y (1) Xtrain = 6 y (2) 4 ..

.

. . . 50 8000 (4) (5) (6) 8150 samples (/5min) 8050 8100 8200 8250 8300 where each row can be considered as an observation and each column is the values of a single variable. 3

Methods

A multivariable statistical technique based on LVs [ 14, 8 ] is developed to predict the future BGC values, where a LV model is developed using the PCA algorithm. Then, the conditional mean values are estimated from the same distribution, and future BGCs are predicted with LVs and incomplete observations. 3.1

Latent Variable Model

Consider a data matrix X (N M ) that contains N observations (rows) and M variables (columns), then a linear combination of an observation x in dataset X can be written as t = p1x1 +: : :+pM xM , where t is a new vector in the same space as x. The fundamental idea behind PCA is to find the loading vector p that maximizes the variation of t, thus the first LV of PCA model can be calculated by solving the following problem: arg max tT t = arg max pT XT Xp

kpk=1 kpk=1 where p is the vector of regression coefficients. Accordingly, the dataset X can be expressed as X = tpT + E with E denoting the residual matrix. We can get more components by solving the following problem:

Traditionally, the successive progress is evaluated by arg max X kpk=1

tpT 2 kXk2

kEk2 100% kXk2 which is referred to as the percentage of explained variation of t and it is fixed as 95% in this study. If the first A (A is much smaller than the rank of X) significant components can summarize sufficiently well the dataset X, then a LV model developed using PCA algorithm can be expressed as

X = t1p1T + t2p2T + ::: + tApTA + E = T1:AP1T:A + E (7) where T1:A = [t1; : : : ; tA] (N A) and P1:A = [p1; : : : ; pA] (M A) are now matrices containing A score vectors and A loading vectors, respectively. For a new test observation z that was not used in the model development but is drawn from the same distribution as observations in X, the scores of the first A significant components of the new observation can then be calculated as 1:A = P1T:Az (8)

Consider a situation where only part of the object z is observed, it is nature to assume the first R variables are measured and the remaining (M R) variables are unobserved, then without loss of generality, we have z = z# z X =

X#

X P =

P # P where z# is the observed variables and z represents the missing measurements. This induces the following partition in X: where X# contains the first R columns of X and X is made up of the remaining (M R) columns. Correspondingly, the loading matrix P can be partitioned as where P # is the submatrix comprising the first R rows of P and the remaining (M R) rows are defined as matrix P . Then, the score vector 1:A of the new observation z can be rewrite as

^1:A = P1#:AT z# + P1:TAz

For a given data matrix X and a new observation z with z# denoting the measured variables that follow the same distribution as observations in X, the missing variables z of z can be estimated as the expected values from the conditional normal distribution: z^ = E z z#; S

Substituting the expression into the estimator of score vector of the new observation yields: ^1:A = 1:AP1#:AT S##

z# S## S #

S# where S = XT X =

N 1

S X and is the square diagonal matrix consisted of the eigenvalues = [ 1; : : : ; H ] in descending order and H is the rank of X.

Finally, the unmeasured variables z in z can then be calculated from the estimated score vector along with the loading matrix: is the covariance matrix of z^ = P1:A ^1:A which yields the future predictions given the past observed glucose values. (9) (10) (11) (12) 3.3

Glucose prediction based on LV model

An accurate glucose prediction model could benefit diabetes management and significantly reduce the risk of hypoglycemia. To predict the future 60 minutes glucose values, the preceding one hour of data is assumed available and marked as z#. Then, the glucose prediction model can be developed and the subsequent future hour of BGC values can be predicted online as follows: 1. The batch data set Xtrain is generated for each subject using the training dataset. 2. While a new observation z# is available: (a) Calculate the similarities between the new object z# and observations in submatrix X#, select the most similar N observations from Xtrain to form a new data matrix X. (b) Develop PCA model using data matrix X as stated in (7). (c) Estimate the score vector of the new observation z# as stated in (11). (d) Predict the next 1 hour’s glucose values as stated in (12).

The above process was implemented in MATLAB 2019b and the future one hour’s BGC values (12 samples) can be forecasted by feeding the testing data to the model. The codes are available: https://github.com/xiaoyu1115/BGLP_2020.git. 4

Results

In the clinical study, the training data from Ohio T1DM dataset with subject ID: 540, 544, 552, 567, 584, and 596 were first divided into two parts: the first 90% of data were used to develop PCA model and the model is tested with the remaining 10% of data to determine the number of observations that should be used in building the LV-based glucose prediction models. Using this approach, the computational burden in the modeling progress can be reduced significantly. Feeding both the processed training and testing data into the modeling process described in Section 3, the future glucose concentrations can be predicted with prediction horizons of up to 60 minutes. The root mean square errors (RMSEs) and mean absolute errors (MAEs) of the 6 subjects are summarized in Table 1. The prediction results are also analysed using Clark Error Grid (CEG) [ 5 ] and the percentages of data distributed in Zone A to Zone D are summarized in Table 1 as well.

v uu 1 RM SE = t N

M AE =

N X (y (i) i=1

N N1 X jy (i) i=1

y^ (i))2 y^ (i)j (13) (14) where y is the BGC values measured by CGM sensors, y^ is the predicted BGCs, and N is the total data points in the testing dataset for each subject.

In Table 1, the RMSE varies from 16.66 to 22.76 mg/dL for the prediction horizon of 30 minutes with an average value of 19.37 mg/dL. The RMSE varies from 26.81 to 38.99 mg/dL for the prediction horizon of 60 minutes with an average value of 32.59 mg/dL. It is reasonable to see this increase in RMSE because the unknown disturbances including meals and exercise that will influence the glucose dynamics significantly. The relationships between the prediction horizon and RMSE or MAE in Figure 2 also indicate that the prediction accuracy decreases with increasing prediction horizon. An average of 87.91% data lie in zone A of CEG which indicates the high clinical accuracy of the model with prediction horizon of 30 minutes. The percentage of data in zone B increase dramatically as the prediction horizon increases to 60 minutes.

The forecasting performance of the model is better for subjects 544, 552, and 596 than subjects 540, 552, and 567. One of the reason is there are less noise and sensor failures contained in the dataset provided by subjects 544, 552, and 596. A filter might improve prediction performance for this case. And many gaps are observed in CGM measurements, which also deteriorate the prediction accuracy of the model in which the interpolated CGM data were used to predict the further BGC values.

The predicted BGC values for subject 544 with prediction horizon of 30 minutes and 60 minutes are shown in Figure 3. For the prediction horizon of 30 minutes, the LV-based model can predict the future BGC values with a comparable high accuracy, and most of the hyperglycemia and hypoglycemia events can be forecasted on time. When the prediction horizon is increased to 60 minutes, the prediction accuracy decreases, but is still acceptable. The prediction values can still provide an insight on the variation of BGC which can be used to tune the insulin infusion rate as a reference. 5

Conclusion

In this paper, an online recursively identified LV-based modeling approach is developed to predict the future BGC values with prediction horizons of 30 and 60 minutes. With the pre-processed dataset, the LV-based model is developed to calculate the LVs using a small submatrix of the training dataset when a new observation is available. The future blood glucose values are predicted as a linear combination of estimated scores and loadings. The proposed model is evaluated with the Ohio T1DM dataset and the results demonstrate the effectiveness of the model. Although the measurement noise is weight-averaged in the LV-based model, it still has a significant influence on modeling and prediction progress. Online denoising techniques would be one of the future study directions that might improve the prediction accuracy. Further, integrating other data fields with the personalized physiological models is a potential approach to improve the prediction performance in the future work.

ACKNOWLEDGEMENTS

The work of Xiaoyu Sun was supported by the China Scholarship Council under grant 201906080136. Funds provided by the Hyosung S. R. Cho Endowed Chair at Illinois Institute of Technology to Ali Cinar is gratefully acknowledged.

References

[1] Arthur Bertachi , Lyvia Biagi, Iva´n Contreras, Ningsu Luo, and Josep Veh´ı, ' Prediction of blood glucose levels and nocturnal hypoglycemia using physiological models and artificial neural networks' , in 3rd International Workshop on Knowledge Discovery in Healthcare Data, KDH@ IJCAI-ECAI 2018 , 13 July 2018 , pp. 85 - 90 .

[2] Marc Breton , Anne Farret, Daniela Bruttomesso, Stacey Anderson , Lalo Magni, Stephen Patek, Chiara Dalla Man, Jerome Place, Susan Demartini, and Simone Del Favero , ' Fully integrated artificial pancreas in type 1 diabetes: modular closed-loop glucose control maintains near normoglycemia' , Diabetes , 61 ( 9 ), 2230 - 2237 , ( 2012 ).

[3] Jianwei Chen , Kezhi Li , Pau Herrero , Taiyu Zhu, and Pantelis Georgiou, ' Dilated recurrent neural network for short-time prediction of glucose concentration' , in 3rd International Workshop on Knowledge Discovery in Healthcare Data, KDH@ IJCAI-ECAI 2018 , 13 July 2018 , pp. 69 - 73 .

[4] Ali CinAr , ' Artificial pancreas systems: An introduction to the special issue' , IEEE Control Systems , 38 ( 1 ), 26 - 29 , ( 2018 ).

[5] William L Clarke , ' The original clarke error grid analysis (ega)' , Diabetes technology & therapeutics, 7 ( 5 ), 776 - 779 , ( 2005 ).

[6] Iva´n Contreras, Arthur Bertachi, Lyvia Biagi, Josep Veh´ı, and Silvia Oviedo, ' Using grammatical evolution to generate short-term blood glucose prediction models .', in 3rd International Workshop on Knowledge Discovery in Healthcare Data, KDH@ IJCAI-ECAI 2018 , 13 July 2018 , pp. 91 - 96 , ( 2018 ).

[7] Meriyan Eren-Oruklu , Ali Cinar, Lauretta Quinn, and Donald Smith, ' Estimation of future glucose concentrations with subject-specific recursive linear models' , Diabetes Technology and Therapeutics , 11 ( 4 ), 243 - 253 , ( 2009 ).

[8] Jesus Flores-Cerrillo and John F MacGregor, 'Latent variable mpc for trajectory tracking in batch processes' , Journal of process control , 15 ( 6 ), 651 - 663 , ( 2005 ).

[9] Iman Hajizadeh , Mudassir Rashid, Kamuran Turksoy, Sediqeh Samadi, Jianyuan Feng, Mert Sevil, Nicole Hobbs, Caterina Lazaro, Zacharie Maloney, Elizabeth Littlejohn , et al., ' Incorporating unannounced meals and exercise in adaptive learning of personalized models for multivariable artificial pancreas systems' , Journal of Diabetes Science and Technology , 12 ( 5 ), 953 - 966 , ( 2018 ).

[10] Roman Hovorka , Valentina Canonico, Ludovic J. Chassin , Ulrich Haueter, Massimo Massi-Benedetti, Marco Orsini Federici, Thomas R. Pieber, Helga C. Schaller, Lukas Schaupp, Thomas Vering, and Malgorzata E. Wilinska, ' Nonlinear model predictive control of glucose concentration in subjects with type 1 diabetes' , Physiological Measurement , 25 ( 4 ), 905 - 920 , ( 2004 ).

[11] Kezhi Li , John Daniels, Chengyuan Liu, Pau Herrero-Vinas, and Pantelis Georgiou, 'Convolutional recurrent neural networks for glucose prediction' , IEEE journal of biomedical and health informatics , ( 2019 ).

[12] Cindy Marling and Razvan Bunescu, ' The ohiot1dm dataset for blood glucose level prediction: Update 2020 '.

[13] John Martinsson , Alexander Schliep, Bjo¨ rn Eliasson, Christian Meijner, Simon Persson, and Olof Mogren, ' Automatic blood glucose prediction with confidence using recurrent neural networks' , in 3rd International Workshop on Knowledge Discovery in Healthcare Data, KDH@ IJCAIECAI 2018 , 13 July 2018 , pp. 64 - 68 .

[14] Philip RC Nelson , Paul A Taylor, and John F MacGregor, ' Missing data methods in pca and pls: Score calculations with incomplete observations' , Chemometrics and Intelligent Laboratory Systems , 35 ( 1 ), 45 - 65 , ( 1996 ).

[15] Sediqeh Samadi , Kamuran Turksoy, Iman Hajizadeh, Jianyuan Feng, Mert Sevil, and Ali Cinar, ' Meal detection and carbohydrate estimation using continuous glucose sensor data' , IEEE Journal of Biomedical and Health Informatics , 21 ( 3 ), 619 - 627 , ( 2017 ).

[16] Mert Sevil , Mudassir Rashid, Zacharie Maloney, Iman Hajizadeh, Sediqeh Samadi, Mohammad Reza Askari, Nicole Hobbs, Rachel Brandt, Minsun Park, Laurie Quinn , et al., ' Determining physical activity characteristics from wristband data for use in automated insulin delivery systems' , IEEE Sensors Journal , ( 2020 ).

[17] Qingnan Sun , Marko V Jankovic, Lia Bally , and Stavroula G Mougiakakou, 'Predicting blood glucose with an lstm and bi-lstm based deep neural network' , in 2018 14th Symposium on Neural Networks and Applications (NEUREL) , pp. 1 - 5 . IEEE, ( 2018 ).

[18] Kamuran Turksoy , Elizabeth Littlejohn, and Ali Cinar, ' Multimodule, multivariable artificial pancreas for patients with type 1 diabetes: regulating glucose concentration under challenging conditions' , IEEE Control Systems Magazine , 38 ( 1 ), 105 - 124 , ( 2018 ).

[19] Kamuran Turksoy , Laurie Quinn, Elizabeth Littlejohn, and Ali Cinar, ' Multivariable adaptive identification and control for artificial pancreas systems' , IEEE Transactions on Biomedical Engineering , 61 ( 3 ), 883 - 891 , ( 2013 ).

[20] C. Undey and A. Cinar , ' Statistical monitoring of multistage, multiphase batch processes' , IEEE Control Systems Magazine , 22 ( 5 ), 40 - 52 , ( 2002 ).

[21] Cenk U ¨ndey, Sinem Ertunc¸, Eric Tatara , Fouad Teymour, and Ali C¸ nar, ' Batch process monitoring and its application to polymerization systems' , in Macromolecular Symposia, volume 206 , pp. 121 - 134 . Wiley Online Library, ( 2004 ).

[22] Malgorzata E Wilinska , Ludovic J Chassin, Helga C Schaller, Lukas Schaupp, Thomas R Pieber, and Roman Hovorka, ' Insulin kinetics in type-1 diabetes: continuous and bolus delivery of rapid acting insulin' , IEEE Transactions on Biomedical Engineering , 52 ( 1 ), 3 - 12 , ( 2005 ).

[23] Xia Yu , Mudassir Rashid, Jianyuan Feng, Nicole Hobbs, Iman Hajizadeh, Sediqeh Samadi, Mert Sevil, Caterina Lazaro, Zacharie Maloney, Elizabeth Littlejohn, Laurie Quinn, and Ali Cinar, ' Online glucose prediction using computationally efficient sparse kernel filtering algorithms in type-1 diabetes' , IEEE Transactions on Control Systems Technology , 1 - 13 , ( 2018 ).

[24] Chunhui Zhao , Eyal Dassau , Lois Jovanovicˇ , Howard C Zisser, Francis J Doyle III , and Dale E Seborg, ' Predicting subcutaneous glucose concentration using a latent-variable-based statistical method for type 1 diabetes mellitus' , Journal of diabetes science and technology , 6 ( 3 ), 617 - 633 , ( 2012 ).

[25] Hong Zhao , Chunhui Zhao , Chengxia Yu , and Eyal Dassau, ' Multiple order model migration and optimal model selection for online glucose prediction in type 1 diabetes' , AIChE Journal , 64 ( 3 ), 822 - 834 , ( 2018 ).

[26] Taiyu Zhu , Kezhi Li , Pau Herrero , Jianwei Chen, and Pantelis Georgiou, ' A deep learning algorithm for personalized blood glucose prediction' , in 3rd International Workshop on Knowledge Discovery in Healthcare Data, KDH@ IJCAI-ECAI 2018 , 13 July 2018 , pp. 74 - 78 .