Development of a Model-Based Noninvasive Glucose Monitoring Device for Non-Insulin Dependent People

Continuous-time glucose monitoring (CGM) effectively improves glucose control, as oppose to infrequent glucose measurements (i.e. using Lancet Meters), by providing frequent blood glucose concentration (BGC) to better associate this variation with changes in behavior. Currently, the most widely used CGM devices rely on a sensor that is inserted invasively under the skin. Because of the invasive nature and also the replacement cost of sensors, the primary users of current CGM devices are insulin dependent people (type 1 and some type 2 diabetics). Most non-insulin dependent diabetics use only lancet glucose measurements. The ultimate goal of this research is the development of CGM technology that overcomes these limitations (i.e. invasive sensors and their cost) in an effort to increase CGM applications among non-insulin dependent people. To meet this objective, this preliminary work has developed a methodology to mathematically infer BGC from measurements of non-invasive input variables which can be thought of as a “virtual” or “soft” sensor approach. In this work virtual sensors are developed and evaluated on 20 subjects using four BGC measurements per day and eight input variables representing meals, activity, stress, and clock time. Up to four weeks of data are collected for each subject. One evaluation consists of 3 days of training and up to 25 days of testing data. The second one consists of one week of training, one week of validation, and 2 weeks of testing data. The third one consists two weeks of training, one week of validation and one week of testing data. Model acceptability is determined on an individual basis based on the fitted correlation to CGM testing data. For 3 day, 1 week, and 2 weeks training studies, 35%, 55% and 65% of the subjects, respectively, met the Acceptability Criteria that we established based on the concept of usefulness. Disciplines Biomechanics | Dietetics and Clinical Nutrition | Endocrinology, Diabetes, and Metabolism | Exercise Science | Kinesiology | Process Control and Systems Comments This article is published as Rollins DK, Beverlin L, Mei Y, Kotz K, Andre D, Vyas N, Welk G, Franke WD. The development of a virtual sensor in glucose monitoring for non-insulin dependent people. Bioinformatics and Diabetes. 2014:1, 19; doi:10.14302;issn.2374-9431.jbd-13-283. Posted with permission. Authors Derrick K K. Rollins, Lucas Beverlin, Yong Mei, Kaylee Kotz, David Andre, Nisarg Vyas, Gregory Welk, and Warren D. Franke This article is available at Iowa State University Digital Repository: https://lib.dr.iastate.edu/kin_pubs/38 Freely Available Online www.openaccesspub.org | JBD CC-license DOI : 10.14302/issn.2374-9431.jbd-13-283 Vol-1 Issue –1 Page No 19 ISSN NO: 2374-9431


Introduction
Recent research suggests that real-time, frequent, glucose monitoring can improve blood glucose control over infrequent monitoring provided through the use of lancet glucose meters for both insulin dependent [1]-[6], and non insulin dependent diabetics [7]. Frequent glucose measurement capability is referred to as continuous-glucose monitoring (CGM); although not really continuous, current devices can deliver on-line glucose measurements as often as every one to five minutes [8]. Nonetheless, this is a substantial improvement over lancet monitoring that only produces a few values (e.g. four values) per day, at best. CGM therefore improves the user's ability to achieve better glucose control by providing frequent, real-time, glucose concentration levels that enables correlation with activity and food consumption. For example, a user is able to see with a high frequency display rate the extent to which the size of a meal affects glucose changes.
Currently, the most widely used and effective CGM devices rely on a sensor that is inserted invasively under the skin. Sensors cost from $35 to $60 and last 3 days to a week. Thus, two significant drawbacks of these devices are comfort and cost [9].
Given these drawbacks, these devices are not widely used except by insulin dependent diabetics that rely heavily on a fast sampling rate for better control. For this reason, these devices are less likely to be used by non-insulin dependent people, including non-diabetic, pre-diabetics and diet-controlled type 2 diabetics.
Hence, the motivation of this work is the development of a useful, non-invasive, subject-specific (personalized), continuous monitoring system in an effort to increase CGM among non-insulin dependent people.
To achieve this goal we seek to develop a low maintenance, high frequency monitoring system with an accuracy that is high enough to be useful for non-insulin dependent people. Moreover, this preliminary work proposes an inferential (i.e., virtual) sensor approach for predicting blood glucose concentration (BGC) from noninvasive inputs. This virtual sensor updates at the same rate as conventional physical sensor CGM devices.
The model is developed from lancet BGC measurements that are obtained at a rate of four measurements per day. Since each sensor is calibrated from user data, the model developed for each person is said to be "subjectspecific." While inferential modeling of BGC has been done by a number of researchers [19]- [24], [29], [31] particularly in type 1 diabetic applications using frequent glucose measurements, this is the first approach that we are aware of that seeks to develop an inferential model for non-insulin dependent subjects using infrequent lancet measurements from the subject's personal lancet glucose meter. Our approach to achieve this goal is to use a novel modeling method to infer glucose concentration using non-invasive input measurements for each subject from variables representing food, activity, clock time [10]- [12], and stress [13], [14].

Methods
The main physical component of this system is a BodyMedia ® armband of the type shown in Fig. 1. This device is a multi-sensor monitoring device that provides accurate estimates of physical activity data using accelerometers, heat related sensors and galvanic skin response (GSR) [15]. GSR is the conductivity of the wearer's skin that varies due to physical and emotional stimuli. For more details see [27], [28]. Given that the armband currently uses complex algorithms (e.g., for pattern recognition) it should also be able to incorporate our proposed BGC prediction algorithm. However, this research is beyond the scope of this article which is focused on the development of the modeling methodology. proposed modeling approach are now described.

The Modeling Approach
The basic objective of this work is the development of a subject-specific "soft sensor" or "virtual" sensor methodology that provides "useful" information to help individuals monitor and control their glucose more effectively than with lancet glucose meters. The most critical and challenging objective in this highly underdetermined problem is that the model must be developed from a BGC sampling rate of only four samples per day. These samples will come from the lancet meter of the subject and the idea is to transform these measurements to a CGM display frequency during the period of the day that the subject is not sleeping.
This virtual sensor approach is an inferential model that is developed from measured variables that are termed inputs. This virtual sensor idea has seen wide applications in process monitoring and control applications in recent years [17], [18] due to advancements in computer hardware, software, and measurement technology. Note that to distinguish the type of sensor, i.e., "virtual" versus "physical," we will use the terms "virtual-sensor" and "physical-sensor." In addition, it should be noted that our use of "monitoring" include both the use of a virtual-sensor or physical- The inputs are shown in Table 1.
The ability to map the available input/output information to accurate sensor measurements depends on the model structure, the model building procedure, and the inferential algorithm that we are calling the

M o d e l i n g S t r u c t u r e
The modeling structure of this application must permit accurate parameter estimation under a small number of sampling times (n), effectively handling several inputs with different dynamic behavior, and mild extrapolation.
The proposed modeling network is what we call the Coupled Dynamic Insulin (CDI) network (its structure is given in Fig. 2). As shown, the first input, meal size (x 1 ), The dynamic functions for G i ,i = 1, …, p, I (the I is for insulin), follow the modeling work of Rollins et al. [16] and are second order differential equations of the form: (Continued on page 23) where   and   are outputs from dynamic blocks G 1 and G I respectively, and   to   are the "coupled" model parameters.
We also use backward difference finite derivative approximation on Eqs. (5) and (6) to give Where  t is the error term and assumed to be independently normally distributed with mean 0 and variance σ 2 for all t, and a i , s are static parameters.
As stated in Rollins et al. [16], the modeling objective is simply to maximize the true but unknown correlation coefficient between measured and fitted BGC. This quantity is represented by  y,ŷ and estimated by r fit . Thus, under this criterion, as a minimum, a model is considered useful, if, and only if, Since the degree of usefulness increases with  y,ŷ , the goal is to obtain the largest (as close to the upper limit of 1) value as possible. Due to the highly complex mapping of the parameters into the response space of r fit , the following indirect criterion is used in obtaining the parameter estimates as described in Rollins et al. [16]. (11) Note that only training data are used to compute SSE under Eq. (11).

M o d e l I d e n t i f i c a t i o n P r o c e d u r e
We use the CDI network with Eqs. (1)-(9) and developed a procedure that can accurately estimate the 3(p+1) dynamic parameters, 4 coupled parameters (  ,  ,   and   ) and the p static parameters even when the number of sampling times (n) is much less than 4p + 7, the total number of parameters. This procedure requires each input to have a separate set of dynamic parameters as uniquely met by the Wiener network but not by other common networks (e.g. such as the Auto Regressive Moving Average with eXogenous (ARMAX) variables network) [16].
Let G f,t = 0 and v i,t = 0 for all i in Eq. (9) except for one value of i = j, i.e., v i = v j ≠ 0, for one value of i, i = 2, …, p. Thus, with only one input variable v i = v j , Eq. (9) becomes a simple linear regression model (SLRM). To distinguish this SLM from Eq. (9), the fitted form is written as (12) Where ŷ i,t = the fitted BGC for the one input i at t; ôi and î are the estimated intercept and slope parameters, re-spectively, for the SLRM for input i ; and ŵ i,t is the SLRM estimate of  i , t .
Note that for fitting the SLRM, only five (5) parameters (the temporary static parameters γ 0 and γ i , and the permanent dynamic parameters τ i , ζ i and τ ai ) are estimated each time which, as necessary, is less than n = 12 for three days of data collection, for example.
In Appendix A, a proof is given to show that for the SLRM, r fit = ry t , v i,t . More specifically, for the SLRM, r fit is determined by only ŵ i,t and not by the static model coefficients, ôi and î Thus, for the SLRM, since v i only depends on the dynamic parameters for input i, one can find the set of dynamic parameters that results in the best r fit for each input i separately (i.e., τ i , ζ i and τ ai .).
We exploit this result by decomposing the modeling problem into separate sub-problems that will be In Step 1, our current procedure is to manually adjust the dynamic parameters one input at a time to find the "best" set of values for each input. Our definition of "best" will be given momentarily. In Step 2, the following reduced form of Eq. (9) is applied: (13) where λ 0 is a temporary parameter only used in Step 2.
This step is the most challenging. With a given set of initial values, either some or all the parameters are estimated simultaneously using an effective nonlinear regression algorithm. This process is the most iterative and time consuming as some parameters are manually

Maximize by Minimizing
Subject to : set and fixed and the rest are estimated using the optimization algorithm. This process is iteratively repeated until no more improvement can be made in r fit .

The only input involved in
Step 2 is food, i.e., meal size.
If an adequate r fit (r fit for training should be positive in agreement with Eq (10)) is not found in this step, the modeling procedure is terminated, and we conclude that the procedure failed to find an adequate model for this subject from the given data sets.
Step The "best" set of modeling parameters is determined for two given scenarios. The first one only uses a Training set of data. In this scenario, the goal is to maximize r fit of the training set. Consequently, the procedure is to reach convergence at the global minimum for the least squares objective criterion. This estimation procedure is "unsupervised" training (note this is a different definition from that of T. Hastie's book [25] and A.J. Izenman's book [30]). The second scenario is when there are both Training and Validation sets of data. In this scenario, the goal is to determine the largest r fit for the Validation data set with a "close" value of r fit for the Training data set. Here, convergence for the Training set may not be reached and the Validation set determines when the iterative process terminates.
Since the Validation results determine when the optimization process terminates, this is a type of "supervised" training. This procedure is used to guard against over fitting, (i.e., fitting BGC behavior in the Training set that is not due to true variation in BGC).
The success of both types of training is evaluated through the use of an additional set of data called the "test set" which had no influence on the model identification process (i.e., the parameter estimates).
The first scenario is used when n is small, say 12, the number after 3 days, whereas the second scenario is used otherwise, e.g., when n is 24, the number after 6 days. Parameter estimation was done using the Excel ® Solver Routine. For the missing data due to removal of armband [32]: if data missing lasts for a short period of time (e.g. no more than one hour of missing data), the missing data were interpolated with the average value of the two sides of the missing data interval. If data missing lasts longer than one hour, we set the missing data to its initial value.

D e v e l o p m e n t o f t h e I n f e r e n t i a l E n g i n e
After obtaining a full set of parameter estimates, the proposed model development procedure has two more refinements. The first one is elimination of any armband al. [16] where only the most recent measurement, at t = t*, is used. This equation, which represents the proposed virtual sensor, is given as: (14) subject to: t > t* and 0 <  < 1, where l is an adjustable constant, y t* is the lancet BGC measurement at t = t*, t = the estimated BGC at time t under the Eq. (9) model, t * = the estimated BGC at time t = t* under the Eq. (9) model, and ŷ t = the virtual (i.e., soft) sensor value for the proposed method at time t. Note that (y t* -t * ) represents that amount of correction and this correction diminishes as time increases based on the value of  which is close to 1. Thus, by the time the next lancet measurement is taken, usually ŷ t ≈ t This means that at t = t*, ŷ t ≈ t * ; at t = t* + Δt, ŷ t ≈ y t* ; and for t = t* + kΔt, with k >> 1 and before the next lancet measurement, ŷ t ≈ t . which had only about 3 weeks of data due to loss data), we have obtained results to support the modeling viability. As just stated, these data sets were collected for another study (see Beverlin et al. [26], [32]).
Modifications had to be made to these data sets for use in this study. First, food quantities, which were in grams of carbohydrates, fats and proteins, had to be converted unavailable, then the nearest value was taken with no more than 4 values used per day. The monitoring period was taken to be from 8 am to 10 pm daily which means that this was the only period that virtual BGC were reported. Thus, the period from 10 pm to 8 am was taken to be a non-monitoring period in order to mimic that monitoring is not required during the sleeping period.
Note that, the original data sets contain meal information in terms of grams of carbohydrates, fats and proteins. The amounts were calculated from self reporting logs of the type and quantities of food eaten.
Hence, the errors of these quantities are likely quite high at times and it is likely that a significant number of meals were not recorded or logged at the proper times.
When we converted the quantities to an index value for meal size for this study (i.e. "1" represents small meal size, "2" for medium size, and "3" for large size), we applied the same conversion equation to all of the subjects. Thus, the quality of food information that we developed our models from in this study is quite poor.
Therefore, since these results are obtained under poor food information they indicate the robustness of the technique to low quality food information.
(Continued on page 28)  As a result, subject 21 and 22 were removed from this research from this point on.

Measures of Performance
Model acceptability will be determined on an individual subject basis given that the models are subject-specific and each individual will only be concerned about model accuracy as it pertains to model developed for them. Thus, this study is evaluated based on the number of subject-models that meet a particular Acceptability Criteria. But we state and justify this criteria momentarily, after we present the statistics they it uses.
The first one is called the averaged error (AE) and is simply the average value of the residuals: (15) where m is the number of terms being averaged.
The second one is called the averaged absolute error and is similar to Eq. (15) except that the absolute difference is used for the term in the summation as follows: (16) A scaled AAE value to adjust for spread is used called the relative AAE (RAAE). This measure of performance is determined by dividing Eq. (16) by the standard deviation of the values used to calculate AAE as follows: (17) RAAE is a relative AAE statistic that accounts for large spread in the glucose variation of subjects. For replicated lancet measurements, the study in Rollins et al. [16] determined RAAE to be about 0.60. Thus, we will assume that a value around 0.60 is comparable to the performance of a glucose lancet meter. However, lancet accuracy or even repeatability can vary widely from individual to individual due to the accuracy of the device, inherent variability in the measurement protocol, and human error. Nonetheless, since this is the only result that likely exists in this type of study (i.e., four weeks of data collection under the protocol of this study) we will use it in our criteria. It is also noted that, given that the models in this study are developed from three discrete levels of food size and not from the three types of consumed quantities, and from a much lower frequency of glucose data when comparing to typical CGMS values (i.e., four values per day versus 12 values per hour), we expected RAAE to be higher and allow for slightly higher values in the Acceptability Criteria.
The last statistic or performance measure is r fit .
Based on the results in Rollins et al. [16] for a type 2 diabetic and in Beverlin et al. [26] for the 20 subjects used here, we set a minimum acceptable value for r fit of 0.40. Using these three measures of performance, the Model Acceptability Criteria (MAC) for this study is given as: and RAAE are 18 mg/dL and 0.75, respectively, the AAE sub-criterion is r fit ≥ 0.5 and the RAAE one is r fit ≥ 0.55. Thus, for this example, the fitted model will meet the MAC if, and only if, r fit ≥ 0.5. As this example illustrate, when r fit is in this range, the MAC is written to require r fit values to be higher than 0.4 with this requirement increasing as AAE and RAAE increase. Note that on the upper boundary of r fit = 0.6, AAE must be < 20 mg/dL or RAAE must be < 0.8, and on the lower boundary of r fit = 0.4, AAE < 16 mg/dL or RAAE < 0.6. The lower boundary was defined based on the results in Rollins, et al. [16] and the upper boundary was established from an examination of fitted models in this work. See Fig. 3 for plots of three fitted models as they compare with CGM measured responses at the limits and middle of the MAC.

Results
The results of this study are given in Tables 3-5. Each   table represents a   Since training stop at convergence under the least squares criterion given by Eq. (11), the remaining days consisted of the test set. In addition, for all these subjects, ζ 1 = 0.2, τ a1 = 0, ζ I = 0.8, and τ aI = 0. This was done to increase the degrees of freedom to estimate the more critical parameter τ 1 , τ I and the four coupled parameters and to simplify the optimization.
The best choice for these values is future research work.
As Table 3 shows, the results indicate that 35% of the cases met the MAC. This is really quite promising as a minimum initial calibration period given that the number of data points, n, used is only 12. In practice, it appears that a significant number of subjects could have successful calibration after three days and as more data are collected this number would grow. This conclusion is supported by the results in the next two tables. Table 4 contains results under Eq. (9) for one week of training, one week of validation and two weeks of testing. As shown, 55% of these cases met the MAC.
In addition, for this group that meets MAC versus all the cases in Table 4  We found that using the armband inputs increases r fit for t * by 0.1 over using just food alone. (These cases are not shown for space considerations). Thus, both food and the armband inputs are to obtain the results presented in this section. The robustness to poorer food quality is supported by similar r fit values for this study as compared to the ones in Beverlin et al. [26] and Beverlin [32] where food quantities were used on these same data sets.

Concluding Remarks
This article presented preliminary work on the   This work applied and improved the methodology from Rollins et al. [16]. It was not the purpose of this work to improve on this approach by investigating the impact of other armband variables, or the time variant nature of model parameters, as well as other model improvement issues. To develop the best model for a specific subject, these issues could be considered but they add to the modeling overhead, which is already very high given the small amount of information. The value of this work lies in that the modeling methodology shows great potential in We have overcome many challenges such as the use of a food index, the lack of initial conditions, frequent and long term removal of the armband and overcome. This includes finding novel ways to improve the accuracy that leads to a higher percent of users meeting the MAC. In addition, the model procedure is quite complex as it requires advanced modeling experience and consists of several steps. One way we plan to improve accuracy is by gaining a better understanding on the bounds of each parameter. To address the model identification issue, we plan to development an estimation algorithm that identifies parameters automatically. This program will reside in the armband and will be used to calibrate the virtual sensor from on-line data. These are areas of future research that we have begun and the results are quite promising.
(A1) Thus, with  i > 0, r fit = r yt,vi,t and for  i < 0, r fit = -r yt,vi ,t This result means that if the correlation of blood glucose concentration (BGC) and v i,t is positive,  i can be set at any positive value and r fit , which will be > 0, will depend only of the behavior of v i,t which is independently controlled by the values of the dynamic parameters associated with v i,t . Conversely, if the correlation of BGC and v i,t is negative,  i can be set at any negative value and r fit will be > 0 and independently controlled by the values of the dynamic parameters associated with v i,t .