Personalized Hybrid Models for Exercise, Meal ... - ACS Publications

Aug 11, 2013 - (6-10, 19) Despite the fact that T1D is prevalent in children and adolescent age groups, majority of the available KDMs are developed u...
0 downloads 0 Views 4MB Size
Download PDF
Article pubs.acs.org/IECR

Personalized Hybrid Models for Exercise, Meal, and Insulin Interventions in Type 1 Diabetic Children and Adolescents Naviyn Prabhu Balakrishnan, Lakshminarayanan Samavedham, and Gade Pandu Rangaiah* Department of Chemical and Biomolecular Engineering, National University of Singapore, Singapore ABSTRACT: Inter- and intrapatient variability in blood glucose (BG) metabolism imposes the need for personalized models in diabetes care. Validation of personalized models using the clinical data of type 1 diabetic (T1D) subjects and inclusion of lifestyle factors like exercise as an input in these personalized models are rarely seen in the literature. In this paper, we have developed personalized BG prediction models with a specialized structure comprising three different classes of models: a mechanistic model for meal absorption dynamics, an empirical model for insulin absorption kinetics and a transfer function model for prediction of personalized BG dynamics. Hence, the proposed model structure is termed a hybrid model (HM). Exercise intensity in these personalized HMs is quantified using a measure called rate of perceived exertion (RPE). The clinical data of 34 T1D subjects are used for model development and two different scenarios of cross validationsame day validation (SDV) and different day validation (DDV). The BG data collected during one of the clinical visits have been used for model development (85−90% of data of this visit) and SDV (the remaining 10−15% data of the same visit); while the data collected during the other day visit have been used for DDV. This is to test the ability of developed HMs in predicting the BG dynamics for a prolonged time period, thereby ensuring their potential to capture intrapatient variability. The fitness and cross validation results of personalized HMs not only show accurate prediction of BG dynamics but also reveal the potential of HMs in capturing both inter- and intrapatient variability.

1. INTRODUCTION Modern health care is gradually abandoning the conventional “one-size-fits-all” therapy and rapidly embracing a personalized approach in the treatment of many alarming diseases. The personalized medicine market in the United States alone is estimated to grow to more than $450 billion by the year 2015.1 Such a major transformation would not have happened without advancements in measurement technology and application of various mathematical modeling techniques that catalyzed an improved understanding of many diseases.2 Although there are many definitions for personalized medicine; the widely accepted definition1 is “the right treatment for the right person at the right time”. Hence, a personalized model that can predict the unique patient dynamics at different periods of time can help physicians in devising safe and efficient treatment strategies in less time and cost ultimately leading to better quality of life for the patients. The two well-established products of personalized medicine that have been prevailing for years in core health care are ABO blood grouping and family history (genetic evaluation tool).3 Incorporating the concept of personalization in diabetes care is important because of the inter- and intrapatient variations in blood glucose (BG) metabolism. Development of personalized BG prediction models for diabetic subjects is practically possible because of the two major things mentioned above the advancements in the BG sensor technology and availability of different BG prediction modeling techniques in the literature. The availability of modern sensors like continuous glucose measurement sensors (CGMS) and GlucoWatch (GW) has made frequent measurements of BG with sampling intervals of 5 min possible. Development of a personalized BG prediction model for each patient can help physicians and healthcare workers: (i) to accurately forecast the BG dynamics © 2013 American Chemical Society

of a patient for various interventions, like insulin, meal, exercise, and stress levels; (ii) to devise optimal insulin, meal, and exercise intervention strategies that can help in maintaining the BG level of patients within therapeutic limits for a prolonged period of time; (iii) to educate the patients and their relatives regarding the adverse effects of protocol deviations on BG dynamics, ultimately leading to increased participation of patients and their relatives in strictly following the physician’s prescriptions, which in turn can prevent various long and shortterm diabetic complications; in fact, such diabetic education is strongly recommended by the International Diabetes Federation as evidence-based counselling.4,5 The mathematical models developed for predicting BG dynamics in diabetic subjects can be broadly categorized into two classes: (i) knowledge-driven models (KDMs) and (ii) data-driven models (DDMs). The former class models (also called as mechanistic or first-principle models) are developed on the basis of the physiological knowledge behind the glucose−insulin regulatory mechanism of diabetics,6−11 whereas the latter class models (also called as empirical or black box or correlation-based models) are developed on the basis of the data only.12−15 Our recent review article gives an extensive summary and analysis of both KDMs and DDMs that have been developed so far to predict the BG dynamics in type 1 diabetics (T1Ds).16 In the case of type 2 diabetics (T2Ds), a review on various KDMs available in the literature can be found in Landersdorfer et al.17 There are also a few studies focusing Received: Revised: Accepted: Published: 13020

February 27, 2013 August 7, 2013 August 10, 2013 August 11, 2013 dx.doi.org/10.1021/ie402531k | Ind. Eng. Chem. Res. 2013, 52, 13020−13033

Industrial & Engineering Chemistry Research

Article

models in clinical settings and daily life has always been doubted in the literature because of the frequent need of technical expertise for estimating network connection weights. Also, both the above-mentioned HMs have neither considered the exercise effects on BG dynamics of T1D children nor employed clinical/real time data in model development and validation. In a recent paper,2 we have proposed a personalized hybrid modeling methodology which involved 4 inputs and 1 output. The input block comprised: (i) a knowledge-driven meal absorption model; (ii) basal insulin flow rates at every 5 min; (iii) an empirically derived bolus insulin kinetic model; and (iv) RPE values at every 5 min. The output is BG readings recorded using CGMS sensors at a sampling interval of 5 min. Relationship between the 4 inputs and 1 output was captured using personalized time series models. We applied this methodology on the clinical/real time data of 12 T1D children and adolescents obtained from one of the exercise studies of Diabetes Research in Children’s Network.34 Even though we proposed a HM methodology which incorporated exercise effects with real data, there were some shortcomings like neglecting the role of personalized factors like body weight (W) in the models of meal and insulin, and improper characterization of the basal and bolus insulin inputs, which are explained in section 3.1 of the present article. The absence of personalization in meal and insulin absorption models has made the earlier proposed HM methodology a partially personalized one. Also, it was only tested on a very small cohort of patients. The major objective of this work is to develop and validate personalized HMs for meal, insulin and exercise interventions using the clinical data of 34 T1D children and adolescents obtained from one of the DirecNet exercise studies.34 Hence, the novelty of this work is not only in introducing the concept of personalization via two different classes of models that incorporate exercise, meal and insulin interventions, but also validating these models using the clinical data of a reasonably large cohort of T1D children and adolescents. We have proposed a modified HM methodology which addresses the drawbacks in our previously proposed methodology.2 The major changes introduced in the revised HM methodology are: (i) incorporation of W into meal and insulin input models; (ii) characterizing the insulin input as the sum of basal and bolus insulin rather than defining as two individual inputs. The identified HMs have been validated for two different scenarios: same day validation (SDV) and different day validation (DDV). The rest of this article is organized as follows. The nature of clinical/real time data used in this work is discussed in section 2 with two subsections for nature of input data (sub-section 2.1) and output data (section 2.2). Section 3 focuses on HM methodology; the results and discussion are presented in section 4. Finally, conclusions of this study are provided in section 5.

on development of DDM to predict the BG dynamics of T2D under free living conditions.18 Because T1Ds are the focus of this work, a brief summary on the state-of-art models developed for T1D is provided here. It is well-known that BG management in T1Ds is only via insulin therapy, as there will be no pancreatic insulin secretion in these patients. Lifestyle factors like meal, exercise, and stress will also play a role in determining the BG dynamics. Most of the KDMs for T1Ds predict variations in BG levels for only two inputs: meal and insulin.6−10,19 Despite the fact that T1D is prevalent in children and adolescent age groups, majority of the available KDMs are developed using adult patient cohorts6,8,10,20,21 with a few exceptions like the model of Chen et al.7,19 and virtual patients of UVa simulator.22 However, both these models have not incorporated the effects of exercise on BG dynamics of children and adolescents. Also, many model-based control studies23−25 generate or use virtual patients from the existing nominal models rather than focusing on development of reliable personalized models. On the other hand, various DDMs based on artificial neural networks (ANNs)12 and time series models26,27 have been developed to predict BG dynamics of T1Ds, for inputs related to meal, insulin, skin impedance, heart rate, energy expenditure, near body temperature, etc., without any emphasis on personalization. Moreover, the variables that were used to quantify exercise in both the classes of models require sophisticated devices or special expensive sensors for their accurate measurement, thus making the use of such models questionable for routine practice. This accentuates the necessity of employing exercise intensity variables that can not only be measured at cheaper cost but also with reasonable accuracy. Rate of perceived exertion (RPE) is a variable that defines the level of exertion an individual feels during exercise and can be self-recorded easily via simple talk or pictorial tests using one of the scales available in the literature.28−30 Despite the RPE’s recognition in exercise settings over the past two decades, it has not been used as a variable for quantifying exercise intensity in most of the existing exercise models related to T1D, except our recent work.2 Because the measurement of RPE does not require the use of any special devices or sensors, consistency of RPE measurements may be questioned. However, clinical studies have proved that the RPE can serve as a reliable and valid exercise marker in children and adolescents, except for those who have cognitive disorders.31,32 Also, subjects and their relatives will have to be provided with appropriate training for using the RPE scales. This will enable the patients to be consistent in their perception of exertion during exercise, which in turn will help in reducing the intrapatient uncertainty when recording the RPE values. Hence, with such proven track record from clinical studies mentioned above and with training from healthcare workers, RPE can be used as a reliable marker in the personalized models of T1D children and adolescents. Apart from the KDMs and DDMs, there is another class of models known as hybrid models (HMs), which is not frequently seen in the BG modeling literature. Mougiakakou et al.14 have developed a hybrid model for four diabetic children by combining the compartmental meal and insulin input models (KDMs) with ANN models (DDMs). The outputs of the input compartmental models served as input to the ANN models, which predict BG dynamics. Recently, Zarkogianni et al.33 employed a similar HM methodology to develop personalized models using the virtual T1D data of UVa simulator.22 However, the practical applicability of such ANN

2. NATURE OF CLINICAL/REAL TIME INPUT−OUTPUT DATA This section explains the nature of real time DirecNet data,34 which have been used in the present study to develop personalized HMs. There are 57 T1D patients in the DirecNet exercise cohort considered for this study. The clinical data of some patients in the DirecNet exercise data set are found to be either not clear or missing. After excluding those patients from the cohort, the data of remaining 34 patients have been selected 13021

dx.doi.org/10.1021/ie402531k | Ind. Eng. Chem. Res. 2013, 52, 13020−13033

Industrial & Engineering Chemistry Research

Article

Figure 1. Distribution of input data in 34 T1D children and adolescents; box plots show maximum, minimum and mean values of each input (half and full filled circles). The red diamonds with blue curves on the left side of each box plot shows the distribution of data among 34 patients. Meaning of x-axis labels can be found within the text description of section 2.1.

Figure 2. Mean ± s.d. values of CGMS recorded BG readings during the two clinical visits. Mean values are represented using red circles, the boundaries of mean ± s.d. are represented as blue lines, whereas the green shaded region indicates the spread of measured BG readings in 34 T1Ds.

for investigation in this paper. Each patient has the data measured during two outpatient clinical visits, with each visit lasting 7−8 h starting from noon (before lunch) till evening (before dinner) with a 75 min exercise period in the late afternoon. The major difference between the two visits is that the basal insulin supply was stopped during exercise in one of the visits (defined as SV in this article), whereas it was continued during the other visit (defined as CV in this article). 2.1. Nature of Input Data. This subsection describes the nature of input data extracted from DirecNet database.34 Exercise, meal and insulin are the three inputs to the models

developed in this work. Figure 1 shows the distribution of input data in the 34 patients. The input data related to meal and insulin can be seen in panels A and B in Figure 1, whereas the exercise data can be found in Figure 1C−F. The W of T1Ds used in this work varies from 39.5 to 92.9 kg with a mean of 62.7 kg; see Figure 1[A] for distribution of W in 34 patients. The distribution of insulin (Ins) to carbohydrate (CHO) ratio for lunch and of the usual prelunch bolus doses of 34 T1Ds are shown in Figure 1B. The amount of CHO consumed during lunch has been calculated based on Ins to CHO ratio and usual bolus dose, as it is mentioned in the 13022

dx.doi.org/10.1021/ie402531k | Ind. Eng. Chem. Res. 2013, 52, 13020−13033

Industrial & Engineering Chemistry Research

Article

Figure 3. Structure of proposed HM with input and output blocks.

DirecNet exercise study34 that the lunch and prelunch bolus during the two clinical visits were given in the same way as that of home settings. In addition to lunch, short meals were given to some patients at different times because of incidences of hypoglycemia before or during exercise. The distribution of variables related to other major meals like breakfast, dinner and bedtime snack are not shown, as we have used only the BG data from lunch to just before dinner in this work. In the DirecNet study,34 the basal and bolus insulin doses were supplied using insulin pumps during both the visits. The total prelunch bolus dose administered to a patient during a particular visit is the sum of usual bolus dose (whose distribution is in Figure 1B) and the calculated correction dose, which has been computed based on two input variables - insulin correction factor (CF) and premeal BG target (PMT). The distribution of both these variables can be seen in Figure 1A. The amount of basal insulin infused for each patient has not been provided directly in the DirecNet study;34 rather they have provided other details like CF, PMT, and total premeal bolus dose injected in a day (obtained from the sum of prebreakfast, prelunch, and predinner doses), In the case of exercise, the 75 min treadmill exercise period comprised 4 phases, with each phase lasting about 15 min with 5 min resting period between each phase. RPEs of all the 4 phases were recorded during both the visit days, denoted as RPE1-SV, RPE2-SV, RPE 3-SV, and RPE4-SV for RPEs

recorded during basal insulin stopped visit day and as RPE1CV, RPE2-CV, RPE 3-CV, and RPE4-CV for RPEs recorded during basal insulin continued visit day. The distribution of RPE (along with distribution among 34 T1Ds) in each exercise phase during both the visit days is shown in Figure 1 C−F. RPE values in 34 patients vary from 0.5 to 9. 2.2. Nature of Output Data. The only available measured output data in DirecNet study is BG, which has been measured using CGMS with a sampling interval of 5 min. The mean ± standard deviation (s.d.) of raw BG data for 34 T1Ds during their two clinical visits can be found in Figure 2. The lower range of measured BG (i.e., mean − s.d.) enter into the zones of mild hypoglycaemia with BG < 70 mg/dl during basal insulin CV (especially during exercise periods), whereas the lower range of BG during basal insulin SV remains within the safer limits.

3. DEVELOPMENT OF PERSONALIZED HYBRID MODELS (HMS) The proposed HM has two major blocks depicting the input and output dynamics as illustrated in Figure 3. The input block has three subsystems that explain the insulin, meal, and exercise dynamics, whereas the output block describes the overall BG dynamics of a patient along with the individual input effect BG dynamics. In the case of input block, the insulin and meal absorption dynamics are explained using reliable empirical and 13023

dx.doi.org/10.1021/ie402531k | Ind. Eng. Chem. Res. 2013, 52, 13020−13033

Industrial & Engineering Chemistry Research

Article

Figure 4. Three main stages involved in development of personalized HMs for 34 T1D children and adolescents.

dynamic data at the sampling interval of 5 min using reliable pretreatment models available in the literature. 3.1.1. Pretreatment of Static Insulin Data. Berger’s insulin kinetic model36,37 has been employed to pretreat the static bolus insulin data (Stage A2 in Figure 4). Usually, diabetic subjects are prescribed with two different types of insulin doses, namely, basal and bolus doses. Both these doses are infused subcutaneously via insulin pumps or syringes. In case of insulin pumps, basal dose will be pumped throughout the day, whereas bolus dose will be injected before meal. In DirecNet exercise study,34 subcutaneous basal and bolus doses were given via insulin pumps. The input data corresponding to basal insulin are readily available in the dynamic form, as the DirecNet T1Ds were infused with basal doses at constant or varying flow rates using insulin pumps. As mentioned in the Introduction of this article, there are some limitations in our preliminary work2 which employed a similar modeling strategy on 12 patients of the DirecNet exercise cohort.34 One important limitation of the methodology in Balakrishnan et al.2 is that the developed HMs were not fully personalized due to the use of generalized meal and insulin absorption models. W of a patient is an important personalization factor in determining the insulin dose prescription and meal absorption, which was overlooked in Balakrishnan et al.2 Hence, in the present study, we have personalized the meal and insulin absorption models using W. Another limitation of the methodology in Balakrishnan et al.2 is not considering the fact that the basal and bolus insulin doses were supplied as a single input using insulin pump in the DirecNet exercise study.34 This

mechanistic models, respectively. On the other hand, the overall and individual effect BG dynamics are explained using personalized transfer function models (TFMs). Eventually, the resulting full model (as illustrated in Figure 3) is defined as hybrid model (HM) because of the involvement of three different classes of models in the same model structure. A more elaborate mathematical description on input and output blocks of the proposed HM can be found in the following subsections. Stepwise methodology for developing personalized HMs is illustrated in Figure 4. The three main stages involved in this methodology are as follows: Stage A, pretreatment of input data; stage B, preprocessing of input-output data; and Stage C, identification and validation of personalized HMs. The substages involved in these three main stages are explained in this subsection. 3.1. Stage A: Pretreatment of Input Data. Insulin, meal, and exercise are the three inputs to the HMs. To develop personalized HMs for BG prediction, it is essential that the input data should be dynamic in nature as the CGMS readings were recorded at the sampling interval of 5 min. DirecNet study34,35 has recorded the input data related to exercise and basal insulin in a form that can be easily represented as dynamic data; however, input details related to meal and bolus insulin are recorded in the form of static variables, namely, Ins:CHO ratio and amount of bolus insulin during different meals (see section 2.1 for details). Unlike BG, it is not easy and cheap to measure the meal absorption dynamics from gut to plasma and subcutaneous insulin absorption kinetics. Hence, the static input data related to meal and bolus insulin are converted to 13024

dx.doi.org/10.1021/ie402531k | Ind. Eng. Chem. Res. 2013, 52, 13020−13033

Industrial & Engineering Chemistry Research

Article

Figure 5. Pretreated input dynamic data of Pt. 32 during basal insulin stopped (subplots [A], [C], and [E] for meal, insulin and exercise, respectively) and continued visits (subplots [B], [D], and [F] for meal, insulin, and exercise, respectively).

unnecessarily increased the number of system inputs (1, 1, 2 inputs for meal, exercise, and insulin, respectively) which led to development of 4 input−1 output time series models, thereby increasing the model complexity. Moreover, many studies in the literature model the insulin input as a sum of basal and bolus doses but not as two independent inputs. Hence, in the present study, 3 input−1 output models using a single input for insulin are developed, as illustrated in stage A2 of Figure 4. The total insulin absorption (UI mu/min) can be represented as UI(t ) =

1 [Ubasal(t ) + Ubolus(t )] W

TDD =

1800 CF

(2)

Next, total amount of premeal bolus dose (TDDbo) injected in a day is calculated from n

TDDbo =

∑ Dbo(m) k=1

(3)

Here, Dbo is the amount of bolus insulin dose given at kth premeal bolus infusion. Finally, daily basal insulin dose (TDDba) is calculated from TDD and TDDbo

(1)

TDDba

Here, Ubasal and Ubolus represent the absorption dynamics of basal and bolus insulin, respectively, and W is the body weight of the patient. As mentioned in section 2.1, the DirecNet exercise study group has not provided the basal insulin data directly in their public access database. Hence, the daily amount of basal insulin infused to each patient is calculated on the basis of the formulas that are usually employed by physicians during insulin prescription. The amount of basal insulin infused in a day has been obtained by substituting the existing DirecNet data like CF, PMT, and total premeal bolus dose injected per day (see Figure 1 for the nature of data distribution in 34 patients) into the usual prescription formulas of physicians. The stepwise calculations are explained below. First, the total daily dose (TDD) of insulin is calculated based on CF by using the relationship that is suitable for Humalog or Novolog formulations

⎧ TDD − TDDbo ; TDD > TDDbo ⎪ = ⎨ TDD ; TDD ≤ TDDbo ⎪ ⎩ 2

(4)

The amount of basal insulin dose at every minute is calculated by dividing TDDba with a conversion factor: (24 h × 60 min). Because the basal insulin is supplied continuously throughout the day, the basal insulin is assumed to be continuously absorbed in the plasma space which can be mathematically denoted as ⎧ t ≥ tba(i) ⎪ D ba (i); Ubasal(t ) = ⎨ ⎪ t < tba(i) ⎩ 0;

(5)

Dba and tba are the vectors of daily basal dose and time, respectively. Bolus dose is infused at a single instant, 5−15 min before each meal. Hence, the absorption of bolus insulin in the plasma space will not be continuous and constant. Eventually, 13025

dx.doi.org/10.1021/ie402531k | Ind. Eng. Chem. Res. 2013, 52, 13020−13033

Industrial & Engineering Chemistry Research

Article

Figure 6. Residual analysis plots of Pt. 34; subplot [A] is the ACF plot indicating autocorrelation of output residuals; subplots [B], [C], and [D] are the CCF plots indicating cross correlation between the inputs (meal, insulin and exercise, respectively) and the output residuals.

respectively. In the present work, personalization is incorporated by using the point estimate values of tmax,G (which has been estimated in one of our recent works on personalized mechanistic models).38 The point estimates of tmax,G obtained in our recent work range from 8 to 150 min for 34 patients. Although there are some limitations in the approach of our recent work on personalized mechanistic models,38 it is better to use the personalized estimates of tmax,G rather than employing the generalized meal absorption model for all patients. The value of f in the meal absorption model is kept at constant value as in Hovorka et al.10 3.1.3. Pretreatment of Exercise Data. Finally, in the case of exercise, there is no need for such pretreatment models, as RPE values were recorded during four phases of exercise (stage A3 in Figure 4). During exercise, RPE is recorded by the patient based on his/her perception of exertion in different phases of exercise. In the nonexercise periods, RPE can be assumed to be 0 as the patients were not involved in any other physical activity. Hence, during rest, the value of RPE is 0; and the RPE values recorded during each phase of exercise are used for every 5 min time interval. The pretreated input profiles of meal, insulin, and exercise of Pt. 32 can be found in Figure 5. 3.2. Stage B: Preprocessing of Input-Output Data. The first step in preprocessing (stage B1 in Figure 4) is removing the means from the original input−output data, denoted mathematically as

the absorption kinetics of the bolus insulin has to be obtained using reliable models in the literature. In this work, the absorption kinetics of bolus dose in plasma space is described using Berger’s kinetic model.36,37 Because Humalog bolus infusion is assumed for all the DirecNet patients considered in this work, the mathematical form corresponding to Humalog alone is presented here. For a more general mathematical representation involving various types of commercial insulin formulations, readers are referred to Chen et al.7,19 n

Ubolus(t ) =

∑

(t − tbo(k))s − 1sT50 sD bo(k) (T50 s + (t − tbo(k))s )2

k=1

× ξ(t − tbo(k)) (6)

In the above equation, Dbo and tbo denote respectively the bolus insulin dose and infusion time of kth premeal bolus (out of n infusions) in a day; ξ(t − tbo) is the unit step function whose value is 1 when t ≥ tbo. Here, s is the absorption rate which varies depending on the type of insulin formulation, and T50 is the time interval to attain 50% absorption of injected insulin given by T50 = aD bo + b

(7)

where the parameters: a and b vary depending on type of formulation given to a particular patient. Values of s, a, and b used in this work correspond to Humalog formulation,19 and they are s = 1.6, a = 5.2, and b = 41. 3.1.2. Pretreatment of Static Meal Data. Mechanistic meal absorption model of Hovorka et al.10 has been used as the meal pretreatment model (stage A1 in Figure 4). This meal absorption model is as follows Umeal(t ) =

1 W

n

∑

UM,sc = UM − UM ; UI,sc = UI − UI ; UE,sc = UE − UE ;

Dm(k)f (t − tm(k))e−(t − tm(k))/ tmax ,G

k=1

× ξ(t − tm(k))

BGreal,sc = BGreal − BGreal ;

(9)

In the above equation, the scaled input data for meal, insulin and exercise are denoted as UM,sc UI,sc and Uex,sc, respectively, whereas the scaled BG data is denoted as BGreal,sc, and the mean values of meal absorption, insulin absorption, RPE values, and real BG are represented as UM , UI, UE , and BGreal , respectively. Next step in the preprocessing is splitting the data set for training and testing (stage B2 in Figure 4), where 85−90% of one day visit data has been used for model identification. Remaining 10−15% of the data is used for SDV, and 100% data of the other day visit is used for DDV. The flowchart in the dashed circle in Figure 4 clearly shows this data split.

2 tmax ,G

(8)

Here, Dm(k) and tm(k) indicate the amount of CHO consumed and intake time of kth meal (out of n uptakes) in a day, respectively; ξ(t − tm) is the unit step function whose value is 1 when t ≥ tm. tmax,G is the time-to-maximum appearance rate of glucose in the accessible glucose compartment (after CHO intake), and f is the fraction of CHO absorbed. Hovorka et al.10 defined tmax,G and f as constants with values of 40 min and 0.8, 13026

dx.doi.org/10.1021/ie402531k | Ind. Eng. Chem. Res. 2013, 52, 13020−13033

Industrial & Engineering Chemistry Research

Article

Figure 7. Measures of Goodness in HM fitting: [A] distribution of calculated msefit values and [B] distribution of calculated %Rfit2 values among 34 patients.

Figure 8. Fitness plots: Profiles of measured CGMS BG readings (red circles), and BGpred (blue continuous line) for 34 patients.

3.3. Stage C: Identification and Validation of Personalized HMs. System Identification toolbox in Matlab has been used to identify the relationship between the preprocessed input and output training data in the form of transfer functions (TF’s). Balakrishnan et al.2 have found that process transfer function models (TFMs) predict the BG values accurately with appropriate gain directions for meal, exercise and insulin when compared to the other classes of linear (like Box−Jenkins and Autoregressive Exogenous input models) and nonlinear (Hammerstein−Weiner and nonlinear ARX models)

time series models. Hence, in the present study, we have opted for TF structure to model the relationship between the pretreated-preprocessed input data and preprocessed BG data. Modeling has been done thoroughly by developing competing HMs which comprise TFs with different number of poles (varying between 1 and 3) and with or without zeros for each input (Stage C1). msefit and %Rfit2 have been computed using the standard definitions of mean squared error and R- squared values, respectively. The competing models are filtered in the subsequent stages (C2−C6) to select the best personalized 13027

dx.doi.org/10.1021/ie402531k | Ind. Eng. Chem. Res. 2013, 52, 13020−13033

Industrial & Engineering Chemistry Research

Article

Figure 9. Number of poles (subplots A−C) and zeros (subplots D−F) in the personalized TFs of meal, insulin, and exercise models.

Figure 10. Personalized parameter values of meal effect TFs for BG prediction.

HM. First, the competing TFMs are filtered by discarding models with relatively high msefit and low %Rfit2 (stage C2). The selected models are then evaluated using residual analysis. If the developed model is reliable, cross correlation function (CCF) and autocorrelation function (ACF) plots should not significantly cross the confidence region. Hence, the models which do not pass the cross correlation (between residuals and inputs) and autocorrelation (within the residuals) tests are discarded in stage [C3]. The ACF and CCF plots corresponding to one of the competing HMs developed for Pt. 34 are shown in Figure 6 for illustrative purpose. It can be observed from this figure that the both the ACF and CCF plots

of residuals are within the confidence region (represented by red dotted lines in Figure 6), indicating that the residuals are not correlated, and hence the competing model of Pt. 34 can be selected and taken to stage [C4] of the modeling methodology. The models which have passed through the residual analysis test are now (stage [C4]) validated for two scenarios (SDV and DDV). The corresponding mse values (mseSDV and mseDDV) and percentage R2 values (%RSDV2 and %RDDV2) have been computed for both the validation scenarios. The validated models are again filtered in stage [C5] based on the mse and percentage R2 values. The successful models are assessed based on complexity in the final stage. The model with relatively 13028

dx.doi.org/10.1021/ie402531k | Ind. Eng. Chem. Res. 2013, 52, 13020−13033

Industrial & Engineering Chemistry Research

Article

Figure 11. Personalized parameter values of insulin effect TFs for BG prediction.

Figure 12. Personalized parameter values of exercise effect TFs for BG prediction.

are developed for all 34 patients. This subsection presents and discusses the results obtained using HM methodology. 4.1. Quality of Fitness. After rigorous modeling stages (shown in Figure 4), personalized HM is identified for each patient. The calculated msefit and %Rfit2 values of all the patients are illustrated as histograms in Figure 7. msefit values for 30 out

fewer parameters (particularly, based on poles and zeros) is selected as the personalized model in the stage [C6].

4. RESULTS AND DISCUSSION The proposed HM methodology in Figure 4 has been implemented on the training data set and personalized HMs 13029

dx.doi.org/10.1021/ie402531k | Ind. Eng. Chem. Res. 2013, 52, 13020−13033

Industrial & Engineering Chemistry Research

Article

of 34 patients are 50% of patients (20 in case of meal and 23 in case of insulin) require zeros. On the other hand, exercise effect models for 18 out of 34 patients do not require any zeros in their TFs whereas TFs with zeros can be observed for 16 patients. Thus, HMs of 34 patients possess TFs of different orders showing the variation of model structure from patient to patient. The resulting personalized TF structure (after system identification) of only one patient is discussed here, for brevity. The personalized values of poles, zeros, gains, and time delays for meal, insulin and exercise effect models of all the patients are illustrated in Figures 10, 11, and 12, respectively. Hence, with all these parameter values, it is easy and simple to write similar equations of TFs for the remaining patients. The personalized TF structure of Pt. 38 involves 3 poles with a zero for meal effect, 2 poles with a zero for insulin effect and 2 poles without any zero for exercise effect. Mathematically, the personalized structure of Pt. 38 is as follows

G E,sc(s) =

(10)

KM(1 + Tz Ms) e(−TdMs) (1 + TPM1s)(1 + TPM2s)(1 + TPM3s) (11)

G I,sc(s) =

KI(1 + Tz Is) e(−TdIs) (1 + TPI1s)(1 + TPI2s)

(13)

Equation 10 represents the combined BG effect model, whereas eqs 11−13 represent the individual effect BG TF’s for meal, insulin, and exercise inputs, respectively. UM,sc, UI,sc, and UE,sc denote the scaled meal, insulin and exercise inputs, respectively, whereas GM,sc(s), GI,sc(s), and GE,sc(s) represent the scaled BG dynamics due to meal, insulin, and exercise inputs, respectively. KM, KI, and KE denote the gains of respective inputs, whereas TdM, TdI, and TdE are the delays corresponding to the three inputs. TPM1, TPM2, and TPM3 are the coefficients of three poles (s = −1/ TPM1, s = −1/ TPM2, and s = −1/ TPM3) and TzM is the coefficient of zero (s = −1/ TzM) of meal effect model. Similarly, the poles and zeros of insulin and exercise effect TF’s can be discussed. Number of poles in a TF determines the order of the model. For instance, in case of Pt. 38, the meal effect TF has 3 poles and hence the patient has a third-order TF model. Similarly, Pt. 38 has second order insulin and exercise effect TF models. The order of each input effect model and parameter values of identified TFs (in Figures 10−12) vary from patient to patient. Hence, in case of HMs, both model structure and parameters tend to vary from patient to patient, indicating the interpatient variability. Values of pole coefficients in all the TFs are positive, which imply negative poles (meaning all poles are in left hand plane, LHP); hence, all the developed TFs are stable, as expected for BG dynamics. Values of poles in all the input effect TFs vary between 0.5 and 99 min except for Pt. 19 and 42, for whom one of the poles corresponding to meal effect TF is >100 min. For Pt.19, the third pole coefficient (TPM3) is >150 min, whereas the second pole coefficient (TPM2) is 125 min for Pt. 42. The personalized gains for meal take a positive value, whereas they take negative values for insulin and exercise, showing that the gain direction is captured correctly in the identified TFs. The time delay in meal, insulin and exercise effect TFs vary between 5 and 100 min, 15−74 min, and 2−75 min, respectively. TFs that have zeros show that only positive zero coefficients (meaning negative zeros) occur in meal effect TFs, while only negative zero coefficients (meaning positive zeros) are found in case of insulin and exercise effect TFs. Positive zeros in insulin effect models occur in 23 patients, whereas they can be seen in exercise effect models for 16 patients. Positive zeros in TFs can result in inverse response which may slow down the system.39 Hence, the insulin and exercise effects on BG can show such an inverse response for patients, whose models have positive zeros. In fact, clinical studies40−42 have reported such an inverse BG response during or at the end of exercise based on the nature and effect of preexercise meal consumption. This could be a reason for the presence of positive zeros as some of the T1Ds in DirecNet exercise study34 were fed with pre-exercise snack meals and/or with snacks in-between the exercise phases based on the measured BG levels. Similarly, the positive zeros obtained for insulin effect TFs (which could lead to inverse BG responses) in some of the patients might be due to the effect of prelunch snacks or meals given to the T1Ds participated in the DirecNet study. 4.3. Cross Validation of the Developed Personalized HMs. The developed personalized HMs for each of 34 patients are cross validated using two scenarios: SDV and DDV. In case of SDV, the personalized models showed a very low mseSDV of 99%), in almost

BGpred,sc = UM,scGM,sc(s) + UI,scG I,sc(s) + UE,scG E,sc(s)

GM,sc(s) =

KE e(−TdEs) (1 + TPE1s)(1 + TPE2s)

(12) 13030

dx.doi.org/10.1021/ie402531k | Ind. Eng. Chem. Res. 2013, 52, 13020−13033

Industrial & Engineering Chemistry Research

Article

the developed HMs possess very good potential for predicting BG data measured immediately after model identification. In the case of DDV too, the HM predicted BG trends are similar to the measured BG trends in each and every patient, as can be seen in Figure 15. There are few minor over and under predictions in some patients, but the direction of BG trend has been captured correctly even in those patients (e.g., see BG trends of Pt. 4, 22, and 23 in Figure 15). The mseDDV values for 27 out of 34 patients are 90%. Recall that the gap between the two clinical visits of DirecNet exercise study varies from 7 to 35 days in different patients. Hence, DDV results indicate that the developed HMs show a promising potential to accurately capture the BG dynamics even after several days of first estimation of the model. This confirms the ability of personalized HMs to capture the intrapatient variations in BG dynamics.

all the patients. Because there is not much variation in both these values, histograms of SDV are not shown here for brevity (only the histograms of DDV are shown in Figure 13). Plots

Figure 13. Measures of goodness in HM validation for DDV: subplots A and B represent distribution of mseDDV and %RDDV2 values, respectively.

comparing the HM predicted BG with that of real data (red circles) for SDV and DDV can be seen in Figures 14 and 15, respectively. The BG profiles of all the 34 patients corresponding to SDV (in Figure 14) show the accurate prediction of real BG values at most of the time points, which is the reason for the very low value of mseSDV. Hence, HMs perform extremely well for the SDV scenario, indicating that

5. CONCLUSIONS Modern healthcare is heading toward a revolutionary concept, referred to as personalized−predictive−preventive−participatory (P4) medicine. This work has incorporated the first two P’s of P4 medicine into the care of T1D children and adolescents, with the help of systems engineering tools like

Figure 14. SDV plots: Profiles of measured CGMS BG readings (red circles), and BGpred (blue continuous line) for 34 patients. 13031

dx.doi.org/10.1021/ie402531k | Ind. Eng. Chem. Res. 2013, 52, 13020−13033

Industrial & Engineering Chemistry Research

Article

Figure 15. DDV plots: Profiles of measured CGMS BG readings (red circles), and BGpred (blue continuous line) for 34 patients.

■

modeling and optimization. A specialized hybrid model (HM) structure involving three different model classes has been proposed to develop personalized BG prediction models using the clinical data of 34 T1D children and adolescents for meal, insulin, and exercise interventions. Rate of perceived exertion (RPE), because of its easy and cheap measurement, has been used as a marker to quantify exercise intensity in these personalized HMs. The developed HMs are validated for two different scenarios depending upon the day on which the data were measured after the personalized model estimation. The modeling and validation results show that the HMs are effective in capturing both intra- and interpatient variability of glucose homeostasis. In our forthcoming works, the predictive capability of the personalized HMs developed in this study will be applied to devise optimal BG control strategies that can prevent the adverse effects of hypo- and hyperglycemia in T1Ds.

■

REFERENCES

(1) The New Science of Personalized Medicine: Translating the Promise into Practice. http://www.pwc.com/mx/es/industrias/ archivo/2011-01-sector-salud-Personalized_Medicine_The_New_ Science_Whitepaper_FINAL.pdf (Last Accessed: 28−01-2013) (2) Balakrishnan, N. P.; Rangaiah, G. P.; Samavedham, L. Personalized Blood Glucose Models for Exercise, Meal and Insulin Interventions in Type 1 Diabetic Children. In 34th Annual International Conference of the IEEE Engineering in Medicine & Biology Society; San Diego, CA, Aug 28−Sept 1, 2012; IEEE Engineering in Medicine & Biology Society (EMBC): Piscataway, NJ, 2012; pp 1250−1253. (3) Teng, K.; Eng, C.; Hess, C. A.; Holt, M. A.; Moran, R. T.; Sharp, R. R.; Traboulsi, E. I. Building an innovative model for personalized healthcare. Cleveland Clin. J. Med. 2012, 79 (Suppl 1), S1−S9. (4) Position Statement−Diabetes Education; http://www.idf.org/ position-statement-diabetes-education (last accessed: 28-01-2013). (5) IDF The Urgent Need: Prevention and Management; http://www. diabetesatlas.org/content/urgent-need-prevention-and-management (last accessed: 26-08-10). (6) Bergman, R.; Ider, Y.; Bowden, C.; Cobelli, C. Quantitative estimation of insulin sensitivity. Am. J. Physiol. Endocrinol. Metab. 1979, 236 (6), E667−677. (7) Chen, C.-L.; Tsai, H.-W.; Wong, S.-S. Modeling the physiological glucose-insulin dynamic system on diabetics. J. Theor. Biol. 2010, 265 (3), 314−322. (8) Dalla Man, C.; Rizza, R. A.; Cobelli, C. Meal Simulation Model of the Glucose-Insulin System. IEEE Trans. Biomed. Eng. 2007, 54 (10), 1740−1749. (9) Fabietti, P.; Canonico, V.; Federici, M.; Benedetti, M.; Sarti, E. Control oriented model of insulin and glucose dynamics in type 1 diabetics. Med. Biol. Eng. Comput. 2006, 44 (1), 69−78.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Fax: (+65) 6779 1936. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank DirecNet for publicly sharing the real time data sets and their timely replies to our queries. The first author thanks the National University of Singapore (NUS) for the financial support for his doctoral studies. 13032

dx.doi.org/10.1021/ie402531k | Ind. Eng. Chem. Res. 2013, 52, 13020−13033

Industrial & Engineering Chemistry Research

Article

(28) Borg, G. A. Psychophysical bases of perceived exertion. Med. Sci. Sports Exercise 1982, 14 (5), 377−81. (29) Robertson, R. J.; Goss, F. L.; Boer, N. F.; Peoples, J. A.; Foreman, A. J.; Dabayebeh, I. M.; Millich, N. B.; Balasekaran, G.; Riechman, S. E.; Gallagher, J. D.; Thompkins, T. Children’s OMNI scale of perceived exertion: mixed gender and race validation. Med. Sci. Sports Exercise 2000, 32 (2), 452−8. (30) Williams, J. G.; Eston, R.; Furlong, B. CERT: a perceived exertion scale for young children. Percept. Mot. Skills 1994, 79 (3 Pt 2), 1451−8. (31) Roemmich, J. N.; Barkley, J. E.; Epstein, L. H.; Lobarinas, C. L.; White, T. M.; Foster, J. H. Validity of PCERT and OMNI walk/run ratings of perceived exertion. Med. Sci. Sports Exercise 2006, 38 (5), 1014−9. (32) Utter, A. C.; Robertson, R. J.; Nieman, D. C.; Kang, J. Children’s OMNI Scale of Perceived Exertion: walking/running evaluation. Med. Sci. Sports Exercise 2002, 34 (1), 139−144. (33) Zarkogianni, K.; Vazeou, A.; Mougiakakou, S. G.; Prountzou, A.; Nikita, K. S. An Insulin Infusion Advisory System Based on Autotuning Nonlinear Model-Predictive Control. IEEE Trans. Biomed. Eng. 2011, 58 (9), 2467−2477. (34) Diabetes Research in Children Network (DirecNet) public datasets: The Effect of Basal Insulin During Exercise on the Development of Hypoglycemia in Children with Type 1 Diabetes; http://direcnet.jaeb. org/Studies.aspx?RecID=161 (last accessed: 28-01-2013). (35) Tsalikian, E.; Kollman, C.; Tamborlane, W. B.; Beck, R. W.; Fiallo-Scharer, R.; Fox, L.; Janz, K. F.; Ruedy, K. J.; Wilson, D.; Xing, D.; Weinzimer, S. A. Prevention of hypoglycemia during exercise in children with type 1 diabetes by suspending basal insulin. Diabetes Care 2006, 29 (10), 2200−4. (36) Berger, M.; Rodbard, D. Computer simulation of plasma insulin and glucose dynamics after subcutaneous insulin injection. Diabetes Care 1989, 12 (10), 725−736. (37) Nucci, G.; Cobelli, C. Models of subcutaneous insulin kinetics. A critical review. Computer Methods and Programs in Biomedicine 2000, 62 (3), 249−257. (38) Balakrishnan, N. P. Development of Predictive Models for Diabetics in Routine Life and Emergency Situations. PhD Thesis, submitted to Department of Chemical and Biomolecular Engineering, National University of Singapore, Singapore, 2013. (39) Stephanopoulos, G., Chemical Process Control: An Introduction to Theory and Practice; Asoke, K. G., Ed.; Prentice Hall of India: Delhi, India, 2001. (40) Guezennec, C. Y.; Satabin, P.; Duforez, F.; Koziet, J.; Antoine, J. M. The role of type and structure of complex carbohydrates response to physical exercise. Int J Sports Med 1993, 14 (4), 224−31. (41) Stevenson, E. J.; Williams, C.; Mash, L. E.; Phillips, B.; Nute, M. L. Influence of high-carbohydrate mixed meals with different glycemic indexes on substrate utilization during subsequent exercise in women. Am. J. Clin. Nutr. 2006, 84 (2), 354−60. (42) Thomas, D. E.; Brotherhood, J. R.; Miller, J. B. Plasma glucose levels after prolonged strenuous exercise correlate inversely with glycemic response to food consumed before exercise. Int. J. Sport Nutr. 1994, 4 (4), 361−73.

(10) Hovorka, R.; Canonico, V.; Chassin, L. J.; Haueter, U.; MassiBenedetti, M.; Federici, M.; Pieber, T. R.; Schaller, H. C.; Schaupp, L.; Vering, T.; Wilinska, M. E. Nonlinear model predictive control of glucose concentration in subjects with type 1 diabetes. Physiol. Meas. 2004, 25 (4), 905. (11) Parker, R. S.; Doyle, F. J., III; Peppas, N. A. A model-based algorithm for blood glucose control in Type I diabetic patients. IEEE Trans. Biomed. Eng. 1999, 46 (2), 148−157. (12) Ghevondian, N.; Nguyen, H., Modelling of Blood Glucose Profiles Non-Invasively Using a Neural Network Algorithm. In Proceedings of the First Joint Biomedical Engineering Society/Engineering in Medicine and Biology Conference; Atlanta, GA, Oct 13−16, 999; IEEE Engineering in Medicine and Biology Society: Piscataway, NJ, 1999; Vol. 2, p 928. (13) Ghosh, S.; Maka, S. A NARX modeling-based approach for evaluation of insulin sensitivity. Biomed. Signal Process. Control 2009, 4 (1), 49−56. (14) Mougiakakou, S. G.; Prountzou, A.; Iliopoulou, D.; Nikita, K. S.; Vazeou, A.; Bartsocas, C. S. Neural Network based Glucose - Insulin Metabolism Models for Children with Type 1 Diabetes. In 28th Annual International Conference of the IEEE Engineering in Medicine and Biology Society; New York, Aug 30−Sept 3, 2006; IEEE Engineering in Medicine and Biology Society: Piscataway, NJ, 2006; pp 3545−3548. (15) Van Herpe, T.; Espinoza, M.; Pluymers, B.; Wouters, P.; De Smet, F.; Van den Berghe, G.; De Moor, B. Development of a critically ill patient input-output model. In 14th IFAC Symposium on System Identification; Newcastle, Australia, March 29−31; International Federation of Automatic Control: Laxenburg, Austria, 2006; pp 481−486. (16) Balakrishnan, N. P.; Rangaiah, G. P.; Samavedham, L. Review and Analysis of Blood Glucose (BG) Models for Type 1 Diabetic Patients. Ind. Eng. Chem. Res. 2011, 50, 12041−12066. (17) Landersdorfer, C. B.; Jusko, W. J. Pharmacokinetic/ pharmacodynamic modelling in diabetes mellitus. Clin. Pharmacokinet. 2008, 47 (7), 417−448. (18) Rollins, D. K.; Bhandari, N.; Kleinedler, J.; Kotz, K.; Strohbehn, A.; Boland, L.; Murphy, M.; Andre, D.; Vyas, N.; Welk, G.; Franke, W. E. Free-living inferential modeling of blood glucose level using only noninvasive inputs. J. Process Control 2010, 20 (1), 95−107. (19) Chen, C. L.; Tsai, H. W. Model-Based Insulin Therapy Scheduling: A Mixed-Integer Nonlinear Dynamic Optimization Approach. Ind. Eng. Chem. Res. 2009, 48 (18), 8595−8604. (20) Fabietti, P. G.; Canonico, V.; Orsini-Federici, M.; Sarti, E.; Massi-Benedetti, M. Clinical Validation of a New Control-Oriented Model of Insulin and Glucose Dynamics in Subjects with Type 1 Diabetes. Diabetes Technol. Ther. 2007, 9 (4), 327−338. (21) Roy, A.; Parker, R. S. Dynamic modeling of exercise effects on plasma glucose and insulin levels. J. Diabetes Sci. Technol. 2007, 1 (3), 338−347. (22) Kovatchev, B. P.; Breton, M.; Dalla Man, C.; Cobelli, C. In Silico Preclinical Trials: A proof of concept in closed-loop control of Type 1 Diabetes. J. Diabetes Sci. Technol. 2009, 3 (1), 44−55. (23) Farmer, T. G.; Edgar, T. F.; Peppas, N. A. Effectiveness of Intravenous Infusion Algorithms for Glucose Control in Diabetic Patients Using Different Simulation Models. Ind. Eng. Chem. Res. 2009, 48 (9), 4402−4414. (24) Liu, S.-W.; Huang, H.-P.; Lin, C.-H.; Chien, I. L. Fuzzy-LogicBased Supervisor of Insulin Bolus Delivery for Patients with Type 1 Diabetes Mellitus. Ind. Eng. Chem. Res. 2013, 52 (4), 1678−1690. (25) Ramprasad, Y.; Rangaiah, G. P.; Lakshminarayanan, S. Robust PID Controller for Blood Glucose Regulation in Type I Diabetics. Ind. Eng. Chem. Res. 2004, 43 (26), 8257−8268. (26) Bremer, T.; Gough, D. A. Is blood glucose predictable from previous values? A solicitation for data. Diabetes 1999, 48 (3), 445− 451. (27) Reifman, J.; Rajaraman, S.; Gribok, A.; Ward, W. K. Predictive Monitoring for Improved Management of Glucose Levels. J. Diabetes Sci. Technol. 2007, 1 (4), 478−486. 13033

dx.doi.org/10.1021/ie402531k | Ind. Eng. Chem. Res. 2013, 52, 13020−13033