Using Grammatical Evolution to Generate Short-Term Blood Glucose Prediction models Iván Contreras1 , Arthur Bertachi1,2 , Lyvia Biagi1,2 , Silvia Oviedo1 Josep Vehı́1,3 1 Institut d’Informatica i Aplicacions. Universitat de Girona, Spain 2 Federal University of Technology – Paraná (UTFPR), Guarapuava, Brazil 3 Centro de Investigación Biomédica en Red de Diabetes y Enfermedades Metabólicas Asociadas (CIBERDEM), Spain ivancontrerasfd@gmail.com Abstract advances of continuous glucose monitoring sensors (CGM). Thus, the popularization of CGM sensors has led to more ro- Blood glucose levels prediction provides the possi- bust and portable devices which has stimulated the availabil- bility to issue early warnings related to ineffective ity of semi-continuous BG measurements which in turn are or poor treatments. Advance notifications of ad- frequently used as data source for predictive modeling in di- verse glycemic events can provide sufficient time abetes. windows to issue appropriate responses and adjust the therapy. Consequently, patients could avoid hy- It is well-known in clinical practice that is complex to perglycemia and hypoglycemia conditions which achieve a tight glycemic control specially since certain pa- would improve overall health, safety, and the qual- tients exhibit large variations in their BG signals. There are ity of life of insulin dependent patients. This re- plenty of factors that influence the blood glucose dynamics port concerns to the application of a search-based and thereby influence glycemic control response. Some of algorithm to generate models able to capture the the factors strongly affecting the glucose metabolism are the dynamics of the blood glucose at a personalized pa- exercise or physical activity, weather conditions, dietary dis- tient level. The grammar-based feature generation turbances, physical conditions, psychological status of pa- allows to build complex empirical models using the tients ([Brusko et al., 2005; Fuchsjäger-Mayrl et al., 2002; data gathered by a sensor augmented therapy, a fit- Mianowska et al., 2011]) together endogenous processes, ness band and a basic knowledge of T1D dynamics. such as circadian cycles [Hinshaw et al., 2013], menstrual Final model solutions provide blood glucose levels periods and pregnancies in women ([Evers et al., 2004; estimations using prediction horizons of 30, 60 and Cramer, 1942]) and other diseases. These varied factors are 90 minutes. often complex to identify, and therefore the prediction of BG values using personalized models is specially important in these scenarios [I. et al., 2016]. Customized models can cap- 1 Introduction ture lifestyle factors which influence the physiologic response The human body requires that blood glucose (BG) levels are of a patient to its carbohydrate intake and insulin dosage. maintained in a narrow range, approximately in the range of Thus, the wide range of variability in the glucose dynamics 70 to 110 mg/dl. BG levels are affected by a large number of of T1D patients makes the generation of predictive models a exogenous factors and, therefore, the pancreas is required to challenging and crucial task. regulate these levels by releasing the insulin and glucagon On the one hand, the treatment of diabetes is conditioned hormones that are secreted by β-cells and α-cells, respec- by a high inter-patient variability which leads to a lack of gen- tively. Type 1 diabetes (T1D) is the consequence of an au- eral models to respond to the particularities of patients. On toimmune attack on β-cells that significantly impairs insulin the other hand, intra-patient variability makes it complex to production. Thus, individuals with T1D fully rely on external generalize models for the glucose response of a singular pa- insulin to manage their BG. tient. The variability points at personalized and dynamic glu- The increasing interest in the improvement of the manage- cose models as one of the best options to implements features ment of this disease and its comorbidities is accompanied by to deal with the treatments variability. At present, intelligent several research efforts focused on therapeutic solutions for algorithms are obtaining a substantial success applying data T1D. One of the most challenging efforts is placed in the ar- driven methods to support advanced analytics and providing tificial pancreas (AP) field. AP refers to an automated system individualized medical aid to patients suffering with diabetes. that combine a glucose sensor, a closed-loop control algo- The incremental repositories of data together with the im- rithm, and an insulin infusion device which are all engaged proved performance of intelligent methodologies to handle together to manage BG and reduce T1D adverse events. AP and process this information have led to the development of has promoted the emergence of increasing research in predic- tools and applications that enhance the effective management tion engines [Cobelli et al., 2011] and its role. Additionally, of diabetes [Contreras and Vehi, 2018]. This report propose it has boosted the commercialization and recent technological the implementation of customizable models for patients using Figure 1: Schematic representation of the method implemented to generate prediction models for blood glucose values an evolutionary computation approach. The article focuses on is feasible to evolve useful models that consider BG read- the critical problem of anticipating BG levels in a short-term ings, meals, and insulin dose information to model BG val- (30 to 90 min). The proposal involves a prediction tool based ues. Later, authors extended the findings by including three on the grammatical evolution method which introduces multi- additional virtual patients and using the root mean squared ple features with the aim of dealing with unforeseen changes. error (RMSE) as the fitness function. The authors tested the clinical significance of the results with an error grid analysis 2 Related Work (EGA) by means of Clarke error grid (CEG) and Parkes er- ror grid (PEG). Other previous studies tested the feasibility BG prediction models can be classified into three different of GE prediction systems based on time series of BG levels subsets: physiological models, data-driven models, and hy- [Contreras and Vehi, 2016; Contreras et al., 2017]. These brid models. First, physiological models are usually gener- studies extended the fore-mentioned research to investigate ated by the experts with wide knowledge and comprehen- the utility of a novel and complementary approach by using sion of insulin, glucose metabolism and other parameters. symbolic regression through GE to evolve personalized BG Second, data-driven models completely relies on BG mea- predictive models that incorporate physiological models as surements and other data inputs. These type of models are part of the inputs. These models included the glucose absorp- typically based on artificial intelligence techniques such as tion rate and the insulin on board model. genetic algorithms, robust filters, fuzzy logic, case reason- ing, auto-regressive models, reinforcement learning, random 3 Materials and methods forests, support vector regression, and artificial neural net- works models. Finally, an alternative architecture involves a Figure 1 shows a schematic representation of the overall combination of the two previous approaches. These models methodology proposed in this study. Initially, we collect are commonly used in a pre-processing stage, and the prepro- the experimental datasets (A). Here we use the Ohio dataset [Marling and Bunescu, 2018] which consists on information cessed inputs enter a data driven model. These type of mod- els are commonly known as hybrid models and some recent from a CSII-CGM therapy and the data from a fitness tracker approaches were examined in previous studies [Balakrishnan band. Next, data was subjected to a preprocessed stage (B) et al., 2013; Estrada et al., 2010; Zecchin et al., 2014]. We where we perform an exploratory analysis and a data clean- redirect interested readers to a more comprehensive review of ing tasks. Next, we perform a feature engineering phase (C), prediction BG models in [Oviedo et al., 2016]. which encompasses tasks to provide additional value to the Previous studies using grammatical evolution (GE) to es- dataset. The most representative transformations involve the timate BG values include the studies [Hidalgo et al., 2014; implementation of the following physiological models (D): 2017] in which a novel customization of BG models for five • The insulin on board (IOB): the insulin that remains ac- virtual patient using GE was first proposed. The incorpora- tive within the body [Wilinska et al., 2005]: tion of medical knowledge into the grammar led to the im- dC1 (t) = u(t) − Kdia C1 (t) plementation of an expression for glucose that considered the dt dC2 (t) (1) previous BG values, carbohydrate intake, and insulin admin- dt = Kdia (C1 (t) − C2 (t)) istration. This involved exploring four different grammars IOB(t) = C1 (t) + C2 (t) and five fitness functions and evaluated all the grammars and where the compartments C1 and C2 have initial values functions with respect to all the patients in terms of average set as 0, u(t) is the insulin dose, and Kdia is a constant error as a performance metric. The results indicated that it related to the duration of insulin action (hs) set as 0.013. • The glucose absorption rate RA(t) (mg/min): carbohy- Table 1: General parameters of the GE implementation and its oper- drate intake of the patient. [Hovorka et al., 2004]: ators Cin Cbio t e(−t/tmax,G ) Parameters Value Parameters Value RA(t) = (2) Population size 200 Tournament size 2 t2max,G Generations 500 Max. Wraps 0 where Cin is the amount of carbohydrates digested, Cbio Crossover prob. 0.90 Codon length 256 is the carbohydrate bioavailability, and tmax,G (min) de- Mutation prob. 0.03 Chromosome length 100 notes the time of the maximum appearance rate of glu- cose in the glucose compartment. • The activity on board (AOB): model based on the total character “|”. Each of the alternatives, known as productions, steps of an individual [Ozaslan et al., 2017]. are composed of a sequence of terminals and non-terminals. These definitions indicate that a non-terminal can be substi- AOB(t) = steps(t)e(−ks t) (3) tuted for any of the productions listed. The quality of the generated solutions depends directly on this structure. The where steps(t) is the total number of steps performed at context free grammar proposed here combines insulin, carbo- time instant t and ks is a constant related to the duration hydrates, BG values and physical activity. Furthermore, we of the effects of physical activity set as 0.0115. have also considered the circadian rhythm of T1D patients After pre-processing stages the dataset will provide crucial and the reliance of generated models from the time. A gram- information to the system training and subsequent validation mar is represented by a 4-Tuple N, T, P, S, N being the non- of the method (D). The system requires the definition of a terminal set, T is the terminal set, P the Production rules for problem specific function (E), which evaluates the solutions, the assignment of elements on N and T. And finally, a start and a customized grammar (F), which defines the structure symbol S which should appear in N. Next we present an ex- of the generated solutions. Solutions will be iteratively com- cerpt of the defined grammar: bined to create new and better solutions aiming to incremen- [Body] → ([Term][Op][Body]) | [Term][Op][Body] | tally improve the quality of solutions and reach a final so- ([Term]) | [Term] lution that minimizes the fitness function satisfactorily (G). [Term] → getG([PrevIni], [Cte], [Op], [Preop], [Exp]) | Once a final solution is generated (H), the prediction model is getRa([PrevIni], [Cte], [Op], [Preop], [Exp]) | evaluated using data from the remaining data base informa- getIOB([PrevIni], [Cte][Op], [Preop], [Exp]) | tion. Next, we describe the complete methodology setup in getAOB([PrevIni], [Cte][Op], [Preop], [Exp]) | detail. getCircadian([OpB], [Cte], [Cte], [Cte]) | [Preop] → sqrt | sin | log | pow | exp | [Preop][Preop] | λ 3.1 Grammatical Evolution Approach [Op] → [OpA] | [OpB] The grammatical evolution (GE) method is a type of search- [PrevIni] → 0 | 1 | 2 | 3 | 4 | 5 | 6 |7 | 8 | 9 [OpA] → + | − based algorithm designed to evolve computer programs or ex- [OpB] → / | ∗ pressions defined by a context free grammar, usually defined in Backus normal form (BNF notation). Using the context Where λ defines the empty set which does not contain any free grammar, GE performs a genotype-phenotype mapping terminals. This grammar is designed to constraint the search process which decodes bit strings to generate programs in space of solutions by using functions that operate the histor- an arbitrary language. GE methodology involves two main ical values of the input signals (insulin, carbohydrates, and design phases. On the one hand, the generation of the con- BG). Additionally, they can be biased by a sinusoidal function text free grammar which is in charge of defining the effective to account for the circadian variations of the patients’ physiol- search space of the problem. On the other hand, the search ogy. Summarizing, the generated solutions are combinations process of the solution, that is, the expression derived from of five expressions, namely ([G], [Ra], [IOB], [AOB], and the grammar which states what should be done. Furthermore, [Circadian]). the separate approach for the search and solution spaces can As other evolutionary techniques, the goodness of the fit lead to the generation of complex phenotypes, as all genetic achieved by the generated prediction models not only relies operators are applied to the genotype. The GE method thus in a grammar definition, but also requires the definition of an becomes an attractive method thanks to its flexibility, closely objective. Through this objective, GE measures the value of associated with the high degree of modularity provided by a each solution of the population and guides the methodology well-structured grammar. towards a final solution. The definition of this objective is The context-free evolutionary grammars is a set of deriva- managed by the so-called fitness function which is the other tion rules expressed in the form: key factor of the GE methodologies. In this paper we use the root mean square error (RMSE), see Equation (4), with the aim to train the predictor according to the aim of the present [non-terminal] → {production1 | ... | productionN } challenge. However, there are multiple possibilities to define Each rule is composed by two key elements, a non-terminal a fitness function in this approach, including fitness functions at the left-hand-side {symbol}, and a definition of the non- considering the clinical assessment of the predictions. We terminal at the right-hand-side {productions}. Each defini- have used two of these metrics to assess the results, the glu- tion involves one or more alternatives that are split by the cose specific RMSE (gRMSE) proposed by [Del Favero et al., Table 2: Mean values of the RMSE, gRMSE and Clarke error grid zones for the 6 patients (PH=30, 60 and 90 minutes) PatientP H gRMSE RMSE CEGA CEGB CEGC CEGD CEGE P (559)30 25.11 20.98 88.6 10.5 0.0 0.8 0.0 P (563)30 22.60 19.36 91.6 7.7 0.0 0.6 0.0 P (570)30 23.92 19.55 93.2 6.4 0.0 0.3 0.0 P (575)30 30.10 24.49 82.4 14.3 0.0 3.2 0.0 P (588)30 23.18 20.45 89.5 10.4 0.0 0.1 0.0 P (591)30 24.08 22.28 76.9 19.6 0.0 3.4 0.0 Average30 24.83 21.19 87.1 11.5 0.0 1.4 0.0 P (559)60 35.19 32.47 74.6 23.0 0.3 2.1 0.0 P (563)60 28.43 27.52 75.5 23.3 0.1 1.1 0.0 P (570)60 29.17 26.45 83.2 15.6 0.1 1.1 0.0 P (575)60 32.77 35.29 63.0 32.5 0.5 3.8 0.0 P (588)60 34.78 31.53 74.6 24.3 0.1 0.8 0.1 P (591)60 33.97 34.77 68.0 29.3 0.2 2.3 0.1 Average60 32.39 31.34 73.1 24.7 0.2 1.9 0.0 P (559)90 45.90 39.81 64.7 31.3 0.6 3.4 0.0 P (563)90 37.49 34.22 71.8 26.6 0.0 1.6 0.0 P (570)90 34.81 31.27 81.3 17.9 0.0 0.8 0.0 P (575)90 41.32 39.78 59.4 33.7 0.3 6.6 0.0 P (588)90 40.06 36.85 71.8 26.3 1.1 0.8 0.0 P (591)90 41.76 35.63 58.5 35.7 0.3 5.4 0.0 Average90 40.23 36.26 67.9 28.6 0.4 3.1 0.0 2012] and the Clarke error grid [Clarke et al., 1987]. a fitted representation of hyper and hypoglycemic events and v an satisfactory RMSE values. Regarding the Clarke error u u1 X N grid, the 30, 60 and 90 minutes GE approaches achieved that RM SE = t ˆ (BG(t) − BG(t)) 2 (4) more than 98%, 97% and 96% respectively of the prediction N t=1 results fell inside regions A and B for the test data, which im- plies that the prediction was safe from a therapeutic point of The operators implemented in this GE approach are the view. Predictive models for patients 559, 575 and 591 obtain modulo operator as a mapping operator, the classic crossover the worst results in all the approaches, while models for pa- by a single point, the integer flip mutation, the selection by tients 563, 570 and 588 achieve the best results in that order. tournament and the elitism. Candidate solutions to a the prob- Future work involves further improving the performance and lem were randomly initialized. The related parameters are safety of predictions by generating sets of models based on defined in Table 1. the determination of individual specific dynamics, lifestyle, and other factors. 4 Results Table 2 summarizes the results of all patients and approaches Advantages of GE over other techniques such as linear involved in this study. The results include values for gRMSE, regression, neural networks or support vector machines are RMSE and percentages by zone from the Clarke error grid. In the flexibility, the modularity and the extensive exploratory addition, Figure ?? shows the predictions performed during power of the method which provides ample room for im- the first day of the testing data for patients 563 and 575. provement. The integration of tools with the potential to pro- vide timely warning of poor or ineffective insulin treatment, which could lead to adverse glycaemic events, is of major 5 Discussion And Conclusions relevance for open and closed loop applications. Continuous This paper has described a methodology based on grammati- short-term predictions of BG levels are plausible, but chal- cal evolution (GE) to generate individualized and customized lenging due to variability, delays in insulin and food absorp- models for BG dynamics. The approach method relies on the tion, and delayed BG measurements. Methods that appropri- generation of a set of rules that determine the search space of ately address such delays and variabilities are likely to pro- solutions for a search-based algorithm. We have implemented vide accurate forecasts which would promote numerous po- a hybrid model using GE and insulin on board, activity on tential applications of benefit to T1D patients. These systems board and glucose rate of absorption models to predict BG would involve features such as, warning alerts of immediately values 30, 60 and, 90 minutes ahead. The experimental re- impending problems, automatic recommendations to prevent sults with PH=30 minutes of are promising since they reflects BG excursions or fail detection in pumps or sensors. Figure 2: Comparison between CGM readings and predictions (PH=30) performed during the first 24-h of the test data of patients 563 (patient with the best average RMSE values) and 575 (patient with the worst average RMSE values). Funding XIV Mediterranean Conference on Medical and Biologi- cal Engineering and Computing 2016, pages 1137–1143. This work has been partially funded by the Spanish Govern- Springer, 2016. ment through contracts DPI2016-78831-C2-2-R, ES-2014- 068289 and by the National Counsel of Technological and [Contreras and Vehi, 2018] I. Contreras and J. Vehi. Ar- Scientific Development, CNPq - Brazil (202050/2015-7 and tificial intelligence for diabetes management and deci- 207688/2014-1). sion support: Literature review. J Med Internet Res, 20(5):e10775, May 2018. References [Contreras et al., 2017] I. Contreras, S. Oviedo, M. Vet- toretti, R. Visentin, and J. Vehi. Personalized blood [Balakrishnan et al., 2013] NP. Balakrishnan, L. Samaved- glucose prediction: A hybrid approach using grammat- ham, and GP. Rangaiah. Personalized Hybrid Models for ical evolution and physiological models. PLOS ONE, Exercise, Meal, and Insulin Interventions in Type 1 Dia- 12(11):1–16, 11 2017. betic Children and Adolescents. Industrial & Engineering [Cramer, 1942] HI. Cramer. The influence of menstruation Chemistry Research, 52(36):13020–13033, sep 2013. on carbohydrate tolerance in diabetes mellitus. Can Med [Brusko et al., 2005] TM. Brusko, CH. Wasserfall, MJ. Assoc J, 47(1):51–55, 1942. Clare-Salzler, DA. Schatz, and MA. Atkinson. Functional [Del Favero et al., 2012] S. Del Del Favero, A. Facchinetti, defects and the influence of age on the frequency of cd4+ and C. Cobelli. A glucose-specific metric to assess predic- cd25+ t-cells in type 1 diabetes. Diabetes, 54(5):1407– tors and identify models. IEEE Transactions on Biomedi- 1414, 2005. cal Engineering, 59(5):1281–1290, 2012. [Clarke et al., 1987] WL. Clarke, D. Cox, LA. Gonder- [Estrada et al., 2010] GC. Estrada, H. Kirchsteiger, and Frederick, W. Carter, and SL. Pohl. Evaluating Clinical R. Eric. Innovative Approach for Online Prediction of Accuracy of Systems for Self-Monitoring of Blood Glu- Blood Glucose Profile in Type 1 Diabetes Patients. pages cose. Diabetes Care, 10(5):622–628, 1987. 2015–2020, 2010. [Cobelli et al., 2011] Claudio Cobelli, Eric Renard, and [Evers et al., 2004] I. Evers, H. De Valk, and G. Visser. Risk Boris Kovatchev. Artificial pancreas: Past, present, future. of complications of pregnancy in women with type 1 di- Diabetes, 60(11):2672–2682, 2011. abetes: nationwide prospective study in the netherlands. [Contreras and Vehi, 2016] I. Contreras and J. Vehi. Mid- Bmj, 328(7445):915, 2004. term prediction of blood glucose from continuous glucose [Fuchsjäger-Mayrl et al., 2002] G. Fuchsjäger-Mayrl, sensors, meal information and administered insulin. In J. Pleiner, GF. Wiesinger, A. Sieder, M. Quittan, MJ. Nuhr, C. Francesconi, HP. Seit, M. Francesconi, kinetics in type-I diabetes: continuous and bolus delivery L. Schmetterer, et al. Exercise training improves vascular of rapid acting insulin. IEEE Transactions on Biomedical endothelial function in patients with type 1 diabetes. Engineering, 52(1):3–12, jan 2005. Diabetes care, 25(10):1795–1801, 2002. [Zecchin et al., 2014] C. Zecchin, A. Facchinetti, G. Spara- [Hidalgo et al., 2014] JI. Hidalgo, JM. Colmenar, JL. Risco- cino, and C. Cobelli. Jump neural network for online short- Martin, A. Cuesta-Infante, E. Maqueda, M. Botella, and time prediction of blood glucose from continuous monitor- JA. Rubio. Modeling glycemia in humans by means of ing sensors and meal information. Computer Methods and Grammatical Evolution. Applied Soft Computing, 20:40– Programs in Biomedicine, 113(1):144–152, 2014. 53, jul 2014. [Hidalgo et al., 2017] JI. Hidalgo, JM. Colmenar, G. Kro- nberger, SM. Winkler, O. Garnica, and J. Lanchares. Data based prediction of blood glucose concentrations us- ing evolutionary methods. Journal of Medical Systems, 41(9):142, Aug 2017. [Hinshaw et al., 2013] L. Hinshaw, C. Dalla Man, DK. Nandy, A. Saad, AE. Bharucha, JA. Levine, RA. Rizza, R. Basu, RE. Carter, C. Cobelli, YC. Kudva, and A. Basu. Diurnal pattern of insulin action in type 1 diabetes: impli- cations for a closed-loop system. Diabetes, 62(7):2223–9, jul 2013. [Hovorka et al., 2004] R. Hovorka, V. Canonico, LJ. Chas- sin, U. Haueter, M. Massi-Benedetti, MO. Federici, TR. Pieber, HC. Schaller, L. Schaupp, T. Vering, and ME. Wil- inska. Nonlinear model predictive control of glucose con- centration in subjects with type 1 diabetes. Physiological measurement, 25(4):905–920, aug 2004. [I. et al., 2016] Contreras I., Quirós C., Giménez M., Conget I., and Vehi J. Profiling intra-patient type i diabetes behav- iors. Computer Methods and Programs in Biomedicine, 136:131 – 141, 2016. [Marling and Bunescu, 2018] C. Marling and R. Bunescu. The OhioT1DM dataset for blood glucose level predic- tion. In The 3rd International Workshop on Knowl- edge Discovery in Healthcare Data, Stockholm, Swe- den, July 2018. CEUR proceedings in press, avail- able at http://smarthealth.cs.ohio.edu/bglp/OhioT1DM- dataset-paper.pdf. [Mianowska et al., 2011] B. Mianowska, W. Fendler, A. Szadkowska, A. Baranowska, E. Grzelak-Agaciak, J. Sadon, H. Keenan, and W. Mlynarski. Hba1c levels in schoolchildren with type 1 diabetes are seasonally vari- able and dependent on weather conditions. Diabetologia, 54(4):749–756, 2011. [Oviedo et al., 2016] S. Oviedo, J. Vehı́, R. Calm, and J. Ar- mengol. A review of personalized blood glucose predic- tion strategies for t1dm patients. International Journal for Numerical Methods in Biomedical Engineering, pages e02833–n/a, 2016. [Ozaslan et al., 2017] B. Ozaslan, S. Patek, and M. Breton. Quantifying the effect of antecedent physical activity on prandial glucose control in type 1 diabetes: Defining exer- cise on board. In ATTD 2017, Paris , France, pages A24– A25, Feb 2017. [Wilinska et al., 2005] ME. Wilinska, LJ. Chassin, HC. Schaller, L. Schaupp, TR. Pieber, and R. Hovorka. Insulin