Blood Glucose Prediction for “Artificial Pancreas” System

The aim of modern science in diabetes therapy is to develop a closed-loop system to control blood glucose (BG) ( “ artificial pancreas ” ). Such a system includes glucose monitor, insulin pump and algorithms of their interaction and blood glucose (BG) dynamics analysis. Current work is devoted to mathematic modeling of BG dynamics, development of BG prediction algorithm and its approbation on clinical data Diabetes Research in Children Network (DirecNet). Prediction algorithm is based on sigma model of BG regulation and is supposed to be used for “ artificial pancreas ” system. The algorithm was tested under condition of possible errors in glucometer and/or insulin pump operation. Root-mean-square error of BG prediction is 15.7 mg/dl. Algorithm corrects 97.5% errors.


Introduction
Diabetes mellitus is one of the most widespread and serious diseases. Diabetes mellitus type 1 means that the pancreas loses its ability to produce insulin. The aim of modern science in diabetes therapy is to develop a closed-loop system («Artificial Pancreas»), which can compensate such an insufficiency. «Artificial Pancreas» is a system combining a glucose sensor, a control algorithm and an insulin infusion device [1]. In the case of diabetes mellitus (type 2), pump insulin therapy is not relevant due to metabolism disorder, whereas in the case of insulin-dependent one (type 1), it is the most effective method of diabetes compensation.
Modern systems for blood glucose (BG) monitoring concentration and insulin pumps for subcutaneous insulin injection have certain potential for raising the level of diabetes mellitus therapy. However, any errors in such a system can cause dangerous consequences for the patient's health: hypo-or hyperglycemia. Some errors occur due to biomechanics of the sensor-tissue area (e.g., mechanical pressure in the sensor implantation zone) [2]. In clinical practice, there are two types of errors associated with glucose monitoring: abrupt fluctuations of BG and short-term decrease in the BG sensor sensitivity. Incorrect functioning of BG sensor can be associated with failure in insulin injection rate caused by tubing occlusion.
Moreover insulin and glucose have delays in absorption and action (about 30-100 min for insulin and 5-15 min for glucose) [1]. Thus, insulin should be injected in the way to prevent hyperglycemia.

Modeling of blood glucose dynamics 2.1. Mathematic models of blood glucose and insulin dynamics
In the human body, glucose is the main and the most multifunctional energy source for metabolic processes. Due to utmost importance of stable blood glucose level, there is a complicated regulation system. The scheme of blood glucose and insulin regulation is represented in Figure 1.
Most BG dynamics models tend to describe physiological processes of glucose and insulin regulation. The simplest model describing BG dynamics is the Bergman's minimal model [3]. It is a system of three ordinary differential equations (ODE) (Eq. (1)): and b 6 are rate constants; and b 5 is glucose threshold level (b 5 ).
The advantage of the model (Eq. (1)) is in proper description of BG and insulin dynamics. However, under some circumstances, it behaves incorrectly. The model describes BG peak with returning to the initial BG level g b .
Many models are based on Bergman's minimal model but take into consideration some more parameters [4,5].
Sturis' model [6] includes two negative feedbacks, describing glucose absorption due to insulin action and insulin secretion as a reaction to BG rising. The model is described by the following ODE system (Eq. (2)): where g in is exogenous glucose infusion rate; i p (t) is plasma insulin; i i (t) is interstitial insulin; f 1 (z) is insulin secretion (function of glucose); f 2 (z) is insulin-independent glucose utilization (function of glucose); f 3 (z) is insulin-dependent glucose utilization (function of glucose); f 4 (z) is insulin-dependent glucose utilization (function of insulin); f 5 (z) is glucose production (function of insulin); x 1 , x 2 and x 3 are variables representing delay process between plasma insulin and glucose production; t 1 and t 2 are time constants of plasma and interstitial insulin degradation, respectively; t 3 is delay time between plasma insulin and glucose production; v p and v i are volumes of plasma and interstitial insulin, respectively; v d is volume of glucose space; and e is rate constant for exchange of insulin between the two compartments.
Engelborghs' model [7] is based on Sturis' model (Eq. (1)) and takes into account time delay of insulin secretion and of BG ascent. Insulin injections can be described in the following way (Eq. (3)) (with an assumption that it equals patient's pancreas insulin): where αf 1 gðtÞ is the fraction of insulin delivered by the pancreas of a normal subject into the circulation to maintain blood sugar at its physiological level, τ 1 is pancreatic insulin production time delay and τ 2 is hepatic glucose production delay. In contrast with Sturis' model, Engelborghs' model is sustainable.
Despite significant international experience in modeling of blood glucose and insulin dynamics, there are a number of individual peculiarities of patient that effect BG like his/her physical activity, stresses, etc. Possible ways of taking into account such a variety of factors are in using of indistinct logic and neural networks. As input data for neural networks, information could be used that can have any influence on BG dynamics. The output data are blood glucose concentration or amount of insulin to infuse to normalize BG.
However, processing of neural network algorithms is even more consuming than systems of ordinary differential equations of high orders.

Sigma model
Direct application of mathematic model for BG dynamics is in its predicting capabilities. On its basis, one can make an algorithm for BG prediction that can be used as a part of "artificial pancreas" system. Mathematic models on the basis of high-order differential equations are rather complicated to be used in a processing unit of a portable glucometer, especially when they are complexificated with additional data about numerous factors effecting blood glucose concentration. That is why simple empirical mathematic models [9] for such calculations are more preferable.
Let us consider independently BG ascent and descent processes. The main factor effecting BG increase is food intake. Food can be characterized by amount of carbohydrates and glycemic index. Thus, food intake can be described by three parameters: m f is food mass, υ is relative carbohydrate coefficient and GI is glycemic index. It should be noted that mentioned parameters have to be estimated for each product independently. Carbohydrate mass influences the amplitude of BG increase, whereas glycemic index defines the rate of glucose appearance in the blood stream.
For BG analysis precise track form during BG ascent is not important, only amplitude and assimilation time. The function of accumulating glucose in the blood can be expressed in the following way (Eq. (5)): where GðtÞ is BG value К, G 0 is initial BG value, GI i ¼ GI i =100 is normalized glycemic index, α is normalization coefficient, τ 1 is time of penetration of all glucose into the blood, t i is time of food intake and τ gluc is time delay between food intake and glucose penetration into the blood.
Function of insulin effect on BG can be described as follows (Eq. (6)): where β is normalization coefficient, τ 2 is time of insulin action spike, t j is time of insulin injection, τ ins is delay between insulin injection and its action, n j is the dose of injected insulin and μ is coefficient that refers to the glucose descent due to one unit of injected insulin.
The model is based on logistic function or sigma function. That is why it is call sigma model.
General equation of BG dynamics that takes into consideration food intake and insulin injection can be described as superposition of functions (Eqs. (5) and (6)): Eq. (7) is the simplest variant of BG dynamics model. However, there are other independent factors affecting BG dynamics: -Glucose elimination by kidneys -Glucose accumulation in the liver in the form of glycogen -Glycogenolysis driven by glucagon -Natural residual insulin secretion by pancreas β-cells

-Gluconeogenesis
In current approach, we do not take into account glucose accumulation and glycogenolysis in the human body. These processes do not change BG rapidly and can be evaluated from blood glucose measurements.
As shown in Figure 2, the model (Eq. (7)) describes ascent and descent of BG track quite accurately. It should be noted that model provides better results in the diabetes mellitus type 1 case, because it is characterized by very little or even no insulin secretion. Sigma model appears to be applicable for BG prediction.

BG profile database
For experimental evaluation of developed mathematic model and the prediction algorithm, we need BG tracks with information about factors that affect BG dynamics. Protocols of Diabetes Research in Children Network (DirecNet) are used for that purpose. DirecNet investigates the potential use of glucose monitoring technologies and its impact on the compensation of type 1 diabetes in children [10]. Investigations with completed protocols were directed to estimate accuracy of constant glucose monitor system (CGMS) [11][12][13][14], especially of FreeStyle Navigator sensor [15][16][17] and evaluate physical activities on BG [18,19] and basal insulin injection [20] influence on hypoglycemia in children with type 1 diabetes. Experimental results are represented as a number of tables, reports about patient anamnesis and received data. Data protocols include detailed description of experiments. Tables contain the following information: BG measurements, insulin injections, initial patient data, calibrating BG tracks, anamnesis, meal, physical activity, dawn phenomenon, blood counts (CBC, hormones, biochemical, clinical BG), medicines, stress, etc. Database includes more than 4000 records about food intake, insulin injections and physical activity and more than 2000 glucose daily tracks. This available data can be used for approbation of mathematic models and prediction algorithm.
DirecNet data was structured, analyzed and filtered; daily BG tracks, provided by different glucose meters and clinical blood analysis, were synchronized ( Figure 3).
In Figure 4, there is a BG track with information about food intake and insulin injection. BG growth and reduction accurately correlate with food and insulin intake time and amount.
Clinical analyses were taken as the most precise. Measurement results provided by noninvasive glucometers FreeStyle Flash and One Touch Ultra (measurements were performed once in a 30-60 minutes) match with each other and in the case of no clinical data, their arithmetic mean value was taken as a true one. Supporting track was formed according to these data. Medtronic CGMS and FreeStyle Navigator data (measurements were performed once in a 1-5 minutes) are less accurate and have noticeable value drift; however, due to higher digitalization, they make it possible to reconstruct precise BG track form. Tracks provided by glucose monitors were corrected in accordance with a supporting track. As a result, a set of BG tracks with higher digitalization (1-5 min step) and accuracy was obtained ( Figure 5).
Obtained data were used for approbation mathematic models and prediction algorithm.

Mathematic model approbation
DirecNet databases give enough information for sigma model approbation. However, it is necessary to define a wide range of parameters for other models of BG dynamics.  In particular, Sturis', Engelborghs' and Bennett-Gourley's models require such parameters as quantities of plasma insulin, interstitial insulin and glucose, pancreatic insulin production, insulin-independent glucose utilization, insulin-dependent glucose utilization, the influence of insulin on hepatic glucose production, distribution volumes for insulin in plasma and interstitial fluid.
These parameters were obtained from BG tracks using a neural network. Bennett-Gourley's model has broader set of parameters in comparison with Sturis' and Engelborghs' models. The purpose of the neural network was to reconstruct Bennett-Gourley's model parameters by an input track. Neural network was trained in the following way. Randomly generated set of parameters for Bennett-Gourley's model (from physiological range of values) was used to reconstruct BG track. This track was included in training sample as a stimuli vector and model parameters as a response vector. After training the neural network, it was used to define necessary parameters for models by analyzing track from DirecNet. Neural network processing was implemented by Matlab Neural Net Toolbox using FeedForwardNet, Elman Recurrent Neural Network and embedded neural network. The best efficiency was demonstrated by training algorithm "trainbr" with 14 neurons in the hidden layer and activation functions "logsig" and "purelin" combination. Average error of model parameter determination was 14%. Elman Recurrent Neural Net provided the best accuracy; RMSE was equaled to 27 mg/dl.
It should be noted that track representation by all models appeared to be more precise when neural network was trained by Bennett-Gourley's model. The model parameters are presented in Table 1 and its description is represented in the work of Tolic et al. [21].
We compared results of track modeling by Engelborghs', Bennett-Gourley's and sigma model parameters with BG tracks for T1DM patients from DirecNet databases. Missing parameters of first three models were defined by neural network; sigma model was built according to parameters specified in experiment protocols. Results are represented in Figure 6. Moments of food intake and insulin injections with corresponding doses are also marked on the plot.
As it turned out, Sturis' model is highly oscillating and has significant dips on the track after food intake. Engelborghs', Bennett-Gourley's and sigma models describe BG track form rather accurately; nevertheless, due to the fact that these models show better results in the case of healthy subjects, it can be noticed that BG tends to return to initial value. Moreover, sigma model is more accurate in describing BG track. Model error in the physiological blood glucose concentration range was equal to 15.7 mg/dl. It makes the model applicable for BG prediction algorithm development on its bases. It should be noted that presented RMSE is average for the whole day. The accuracy is higher for short-term prediction (30 min) and model iterative correction.
The disadvantage of the model is its ability to predict only growth or reduction of BG induced by food intake or insulin injection, respectively. As there are a lot of other BG-affecting factors [22], sigma model can be extended and improved but is applicable for prediction algorithm development as is.

Predictive algorithm
BG prediction algorithm is needed for adjusting the rate of insulin injection and detection of insulin pump or glucometer failure.
Algorithm is based on sigma model [23] and calculates BG dynamics for 2 hours in advance. Initial BG level is measured by invasive glucose meter in the fasting state. On receiving new data (from glucose monitor or from patient), a new prediction is made and previous predictions are erased. After each measurement made by BG monitor, the model analyses received data. Incorrect system operation cases and algorithm response on that are described below. If there are no errors in measurement process, new BG dynamic prediction is performed on the base of received data. Linear shift of predicted track toward measured value is calculated as follows (Eq. (8)): where δ gluc is glucose monitor value error and δ m is averaged error of mathematic model for the period between two measurements.
In order to evaluate measured BG adequacy, we consider two types of prediction: short-term and long-term predictions, to detect local errors and systematic ones, respectively.

Error detecting and correction
The first type of errors is related to possible spikes due to BG monitor software or hardware failure during BG measurements. Such kind of errors is revealed to single value of glucometer data that differ significantly from predicted one. Threshold for taking measured value as falling can be expressed as (Eq. (9)) where δ Σ is the total tolerance and δ tol is the theoretical tolerance of BG data deviation from values obtained by the model when accurate measurement and modeling occur.
When such an error takes place, algorithm ignores the data from glucose monitor and considers predicted data as true. Error signalization is not obligatory in this case. During the period of seven measurements, only three falling out values are allowed. If incorrect values appear, too often algorithm makes a report about possible system failure. Algorithm operation example with such type of errors is represented in Figure 7.
The second error type is connected to damage (split or disconnection) or cross clamping of insulin pump infusion set. Because of undelivered insulin, BG systematically increases. If longterm prognosis differs from glucometer data with insulin being infused, the algorithm sends a signal about possible error and suggests patient to check an infusion set ( Figure 8).
It also may happen that a patient for some reason did not enter meal data. In this case, as well as in the previous one, significant deviation of glucometer measurements from long-term prediction occurs. If insulin was not infused, algorithm sends a signal about the error ( Figure 9) and suggests patient to confirm that he/she didn't have a meal during last 60 minutes.
In the case of sharp reduction of sensor sensibility, insulin pump breakdown, or other system operation failure, systematical deviations occur and there is no possibility to eliminate them by prior mentioned methods. Then algorithm reports about system failure, turns off automatic insulin infusion and suggests a patient to consult his/her doctor.

Insulin dosage calculation
The aim of insulin pump is automatic insulin infusion in two regimes: basal and bolus. The pancreas produces insulin in advance so that BG would be in certain range of values. There are several schemes for individual insulin injections and in most cases, it should be injected before BG meal. Rapid-acting insulin should be infused straight after the meal but before blood glucose increase. Moreover, patients often need basal insulin infusion twice a day, but with insulin pump, it is infused continuously [24]. Bolus insulin dose is calculated in accordance to meal composition and product glycemic index.
Using such data and sigma model, one can calculate BG dynamics and evaluate dose and time of insulin injection to make BG normal after food assimilation. As a rule a standard value is used, for instance, 80 mg/dl. As sigma model is sensible only to food intake and insulin injection, BG with t ! ∞ tends to a constant (Eq. (10)): In order to reach desired BG value G norm (Eq. (11)), it is necessary to infuse the following insulin volume n norm (Eq. (12)): Figure 9. Algorithm operation when patient didn't enter meal data.
Blood Glucose Prediction for "Artificial Pancreas" System http://dx.doi.org/10.5772/67142 Insulin infusion time depends on type of insulin used in the insulin pump. Selection of insulin infusion pattern is performed according to definite kind of insulin. The pattern is selected to minimize BG deviation from normal range of values. It is important to prevent hypoglycemic values of BG, because lack of insulin can be compensated by additional bolus, but extra BG reduction caused by insulin over-infusion cannot be balanced automatically.
Predicting track is shifted depending on measured BG values. As a result, it may happen that real BG value after some time will be lower than it was initially predicted. Meanwhile a new prediction will contain a zone where predicting value will appear in hypoglycemia area. Insulin and food have assimilation delays. If insulin of extra value is infused, it will be impossible to register and correct appropriate BG changes and a patient will have hypoglycemia after a while. In order to avoid hypoglycemia, it was offered to infuse only 75% of bolus value at once and the rest after increase of BG.

Concept of predictive algorithm work principles
General concept of predictive algorithm operation is represented by the flowchart (Figure 10).

Experimental tests
Algorithm was approbated using BG tracks from DirecNet database.
Algorithm successfully detected cases of meal data absence. In the case further data input and algorithm operation were impossible to perform. As all tests were performed under the medical supervision, DirecNet protocols have no data about possible insulin infusion system failure.
However, measured BG tracks had many spikes. Results of such kind of error elimination are represented in Table 2.
Algorithm corrected 97.5% of BG monitor errors. Meanwhile only 3.7% of correct BG monitor data was detected as errors. Generally, measurements were taken as errors in the ranges of BG sharp ascent and descent due to food intake or insulin infusion, respectively. That is explained by high track differential. Having that, little variations of delays in food uptake and in insulin infusion lead to significant difference of predicted and real BG values. Uncorrected errors of BG appeared, when spikes had deviated toward the predicted track and then had fallen into algorithm's correct values zone. Along with real errors of BG monitor, the synthetic errors were included in tracks. Approbation of the  algorithm confirmed its sensibility to mentioned situations; however, they rarely happen in practice.

Results and discussion
Sigma model rather accurately describes BG dynamics in patients with diabetes mellitus. For initial model approbation, only food intake and insulin injection data were taken into account.
Other factors affecting BG dynamics can be included in the model after additional research. In this case, more precise clinical data will be needed for mathematic model approbation. Nevertheless, the advantage of sigma model is that it does not require a lot of physiological parameters of patient and can be used for BG prediction as is.
Described model approbation shows that RMSE of Bennett-Gourley model is 18.21 mg/dl. In spite of its simplicity, sigma model has the lowest RMSE among all tested models-15.7 mg/dl. This model contains no differential equations and requires minimal computing power in "artificial pancreas" technical embodiment.
Developed predictive algorithm is based on the sigma model. The algorithm estimates the difference between theoretical and experimental tracks with possibility of further correction. BG prediction algorithm is needed for detection of incorrect measurements and correction spikes on BG track. That allows patient to avoid incorrect insulin infusion or glucose intake and reduce possible risks.
Moreover, the algorithm alerts about closed-loop BG control system failures. It allows patient to detect damage or cross clamping of insulin pump infusion set and actualize meal data. Meanwhile patient at his own risk can continue system operation.
The algorithm allows making glucometer-insulin pump system automatic and its work more physiological. Due to insulin doses and injection time calculation, it is infused in advance (before patient BG ascent). Automatic system avoids errors when patient calculates insulin doses and injects it himself.
The algorithm was tested using 52 tracks with full information about the patient's schedule of food intake, insulin injections, etc. Thus, the algorithm suggested in this work makes it possible to avoid 97.5% cases of incorrect operation of the BG monitor. In this case, the RMSE of the calculated blood glucose concentration values was 11.06 mg/dl.