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

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-sensor although virtual sensors do not measure the process variable being monitored directly. This major challenge in this work is the frequency of BGC data for model building (in this research, 4 times per day) is much less than the virtual measurement rate of 5 minutes. This limitation means that the information available for model identification, i.e., parameter estimation, is quite limited and could thus, severely impact accuracy.

The information for the development of a virtual senor comes from two sources -- the response data set and the input data set. Since the information content of lancet BGC is quite limited, the proposed approach strongly relies on the input data set for information on glucose behavior. More specifically, this data set consists of meal size with three levels, six (6) variables from the BodyMedia armband, and the time of day (TOD) in minutes on the 24 hour clock. The inputs that we selected for this study from the armband are those selected by Rollins et al ¹⁶. We eliminated near body temperature as we determined it was not contributing significantly to glucose behavior for any of the subjects. 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 “Inferential Engine.” The model structure consists of the mathematical functions and the network that tie these functions together. The model building (i.e., identification) procedure is the process of using input/output information to estimate the values of unknown parameters in the mathematical functions. The Inferential Engine is the equation used to obtain the virtual senor measurements at the desired sampling frequency. This equation represents input selection, parameter estimates, and the use of lancet glucose measurements to enhance reliability. The purpose of this section is to describe these three components of the proposed technique in detail.

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 Figure 2). As shown, the first input, meal size (x₁), enters both a linear dynamic food block (G₁) and a linear dynamic unmeasured insulin block (G_I). The output from the unmeasured dynamic insulin block (v_I) enters a pseudo blood insulin block which is coupled with the food block (G₁) which produces the dynamic food input (v₁) to the pseudo BGC block. Then, the unmeasured output from the coupled food block is the dynamic glucose (G_f) input due to food consumption. Each of the other inputs (e.g. inputs 2-8) enters a separate linear dynamic block and the outputs from these blocks are collected into non-observable variables (v_i) and together with G_f are passed through a static block which can be any type of function. The CDI model simulates the process where food digestion is responsible for the rise of blood glucose after each meal, while the secretion of insulin is responsible for the fall of blood glucose level a period of time later after the meal. The CDI network is defined by the attributes of allowing separate dynamic behavior for each input and the use of variables for unmeasured insulin generation (v_I) and unmeasured blood insulin concentration (I). To our knowledge, this is first application of unmeasured pseudo insulin in modeling blood glucose concentration. This idea is a key reason for the success of our modeling approach in this application of infrequent BGC measurements.

Figure 2.A graphical representation of the Coupled Dynamic Insulin (CDI) network.

The dynamic functions for G_i ,i = 1, …, p, I (the I is for insulin), follow the modeling work of Rollins et al. ¹⁶ and are second order differential equations of the form:

(1)

where x_i(t) is the i^th input,i varies from 1 to p, p is the total number of inputs, t_ai is the lead parameter, t_iis the time constant, and z_i is the damping coefficient, with x₁ = meal size input variable, x_i, i= 2, . . ., p-1, are armband input variables, and v_p is the TOD input variable.

Using backward difference finite derivative approximations, Eq. (1) gives (Rollins et al., ¹⁶)

(2)

with

(3)

such that w_2,i. = 1 - d_1,i - d_2,i - w_1,i . This constraint is used to impose a unity gain restriction for the linear dynamic blocks. Δt is the sampling time for the inputs. In the Laplace domain, the linear dynamic functions are

(4)

Note that the number of dynamic parameters associated with each input is three. This small number is a strength that we exploit to obtain parameter estimates under limited sampling, as discussed below. The CDI model for food alone is represented by the following coupled Eqs. (5) and (6):

(5)

(6)

where v₁and v₂ are outputs from dynamic blocks G₁ 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

(7)

(8)

Note there are four additional parameters. (α₁,α₂,α₃ and α₄)that need to be identified.

The function f(V) is called “the static function” and is a function of all of inputs. This function can theoretically be of any form. For effectiveness under mild extrapolation and minimum parameter estimation (as discussed below) we have chosen a first order linear regression model of the form:

(9)

Where _t is the error term and assumed to be independently normally distributedwith mean 0 and variance σ ² for all t, and a_i ^,sare static parameters.

As stated in Rollins et al. ¹⁶, 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,

(10)

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. ¹⁶.

Note that only training data are used to compute SSE under Eq. (11).

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) ¹⁶.

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; ŷ_oi and ŷ_i are the estimated intercept and slope parameters, respectively, for the SLRM for input i; and ŵ_i,tis the SLRM estimate of v_i,t.

Note that for fitting the SLRM, only five (5) parameters (the temporary static parameters γ₀ and γ_i, and the permanent dynamic parameters τ_i, ζ_iand τ_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,tand not by the static model coefficients, _oi andŷ_i 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 bestr_fit for each input i separately (i.e., τ_i, ζ_iand τ_ai.). We exploit this result by decomposing the modeling problem into separate sub-problems that will be identified in 3 steps: 1. the dynamic parameters for each input i, i = 2, …, p, under Eq. (12) (five parameters are estimated for each i); 2. the insulin and food dynamic parameters under Eqs. (7),(8), (13) (eleven parameters are estimated; one temporary intercept parameter, four coupled parameters for initial values to be used in Step 3, and six permanent dynamic parameters); and 3. the permanent static and coupled parameters with all the inputs included under Eq. (9) (at most p + four coupled parameters are estimated).

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 λ₀ 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 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 3 is completed in one estimation trial when the dynamic and couple modeling parameters from Steps 1 and 2 are used since Eq. (9) becomes a first order linear regression model. However, if one desires to estimate the couple modeling parameters also, Eq. (9) becomes a nonlinear regression model and the estimation process can be more challenging. Note that, for example, with n = 12, or three days of data collection and p = 8, at most twelve parameters are estimated in Step 3, which is still not exceeding n. At the end of all three estimation steps, with p = 8, 4p + 7 or 39 total parameters have been uniquely estimated in a highly nonlinear modeling problem from at least n =12 or three days of data collection, for example. This ability is a critical novelty and a powerful benefit of this approach.

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 ²⁵ and A.J. Izenman’s book ³⁰). 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.

Successful model identification relies on effective selection of initial conditions and starting values for model parameters and the dynamic inputs (i.e., the v_i’s). The following procedure is given under a protocol where the armband is worn nearly 24 hours a day and removed only for showering. The initial steady state is chosen during a period of slow change, commonly early in the morning. The set of initial values in our procedure are τ_i = 1.1 ζ_i = 0.9, α₂= 20, α₄= 0.1, and all other initial parameter values are equal to zero. The initial values for the v_i’s have to be determined iteratively. When the dynamic parameters are set to values so are the v_i’s as shown by Eq. (2). Our procedure is to set the initial values of the v_i’s to their average values over the training data. These values are to remain fixed during estimation and changed after estimation of dynamic parameters. The estimation process for a set of dynamic parameters is completed when the “best” r_fit is obtained with initial values of v_i’s close to their average values for the training data.For the missing data due to removal of armband³²: 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.

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 inputs that adversely affect the value of r_fit. This is done by setting each, and only one, a_i (for i = 2, …, 7) to zero at a time, and observing r_fit. If r_fit increases for the Training set in Scenario 1 or for the Validation set for Scenario 2, this input is removed. After this process is completed for each input, all the remaining static parameters are estimated under Step 3 for a final time.

The final refinement involves the use of lancet glucose to help to reduce model bias. Since these measurements are infrequent and are not measured at a constant rate, it is not possible to build a correction model based on the correlation of residuals. The correction equation that we use comes from Rollins et al. ¹⁶ 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*, the estimated BGC at time t under the Eq. (9) model, 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* —}) 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≈ This means that at 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≈ . That is, at the time of the lancet measurement, the proposed virtual monitor would display a value close to , the next value would be close to the lancet measurement, and as time proceeded, the lancet value would have less corrective influence as the predictor would rely more on the model to infer BGC. When two sets are used to estimate model parameters, can be set to give the most accurate values in the validation set. When only a training set is used to estimate the model parameters, a default value can be used based on results from modeling several subjects.

For the proposed method, the development of a virtual-sensor requires 4 lancet measurements per day spread as evenly as possible over the time the subject is awake in about a 14 hour period. We did not have access to data meeting this requirement. However, from a previous study, we had physical-sensor CGM data sets which were collected with Institutional Review Board (IRB) approval, and the data sets were used to develop and evaluate the methodology. Thus, these data sets played two roles. First, for each subject, they played the role of a surrogate person, i.e., the true BGC for the purpose of evaluation. Secondly, they played the role of the lancet sampled data, i.e.,, the data used to build the virtual-sensors.

Using 22 test subjects (see Table 2) with 4 weeks of data collection (in most cases and slightly under 4 weeks in other cases except for Subject 1 and 8 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. ²⁶, ³²). 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 to a food index representing meal sizes with 0 for no meal, 1 (two time stamps) for a small meal, 2 (three time stamps) for a medium meal and 3 (four time stamps) for a large meal. In practice the time stamps will be entered by the user pressing the time stamp button on the armband at the start of a meal. The conversion we used was based on the grams of carbohydrates only with less than 20 grams being a small meal, more than 100 grams being a large meal and all other amounts considered a medium size meal. Secondly, infrequent BGC measurements were not obtained from a lancet meter but converted from a continuous glucose monitoring system (CGMS) at a sampling rate of only four values per day at particular and fixed times (i.e. only 4 values per day out of the continuous readings from CGMS were used) to mimic infrequent lancet sampling. CGMS values were taken only at 8 am, noon, 4 pm and 8 pm; if data were 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.

Before evaluating model acceptability, subject 21 and 22 were rejected due to poor food information. As a result, subject 21 and 22 were removed from this research from this point on.

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. ¹⁶ 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. ¹⁶ for a type 2 diabetic and in Beverlin et al. ²⁶ 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:

(18)

As the MAC shows, a model with an r_fit of at least 0.6 is considered acceptable based on this value alone. However, if this value is between 0.4 and 0.6, the fitted model must meet a certain level of performance on both r_fit and AAE or both r_fit and RAAE. For example, if AAE 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. ¹⁶ and the upper boundary was established from an examination of fitted models in this work. See Figure 3 for plots of three fitted models as they compare with CGM measured responses at the limits and middle of the MAC.

Figure 3.Graphical examples during the testing period for three subjects meeting the MAC: Subject 2 (strongly meeting the MAC); Subject 6 (weakly meeting the MAC) and; Subject 20 (moderately meeting the MAC).

This article presented preliminary work on the development of a virtual sensor for BGC with the objective of developing a noninvasive CGM system that could increase CGM among non-insulin dependent people. This device would require users to wear a readily available armband monitor and manually entering meal sizes through the use of a button on the armband. This device would require four (4) lancet measurements per day as most current invasive CGMSs require.

The modeling methodology presented in this work is quite powerful. It takes on the challenge of modeling BCG in a highly complex, non-linear, multiple-input, highly underdetermined problem. As illustrated in this work, it is able to develop useful multiple-input dynamic models for BGC under free living, outpatient, data collection from just four glucose measurements per day and from as little as three days of data. In addition, these results are achieved with minimal food information of only three discrete levels. This ability stems from a number of innovative ideas to overcome several challenges in this complex modeling problem as follows. First, the use of the coupled structure allows for the inclusion of inferential blood insulin concentration and leads to insulin and glucose interaction in the blood. This structure is a significant advancement over a straight Wiener network and contributes significantly to the accuracy and ability to obtain adequate fitting for acceptable model usefulness. Secondly, the result in Appendix A provided the knowledge that produced the idea to decompose the modeling problem into a dynamic part and a static part. Added to this idea is the inspiration of determining the dynamic parameters for each input, one input at a time. Once the dynamic parameters are determined for each input, they are fixed. Note that, from the use of a validation set we are able to control over-fitting and by controlling r_fit to be about the same in the training set and validation set for each input separately, we have found that this helps the final r_fit in all the data sets (Training, Validation, and Testing) to be quite similar. After obtaining the dynamic parameters, the low number of static parameters is then obtained separately as a linear regression model. Thirdly, as the results show, the correction provided by Eq. (14) contributes strongly to the accuracy of the proposed method in the case of continuous glucose monitoring (CGM). While this correction does contribute significantly to the reduction in bias, it also contributes majorly in the reduction in AAE. This can only occur if there is a significant positive correlation for the fitted response, as the correction brings it close at the times infrequent measurements occur, but the correlation determines the direction from these points. If the correlation is not positive, the trend will be in the wrong direction and accuracy would suffer tremendously. This is why subjects must have a significant positive correlation for acceptability. In practice if this is not achieved, the device would simply give a calibration error and not report measurements. Lastly, we developed the new Model Acceptability Criteria (MAC) which instead of evaluating separate aspects of performance such as correlation or bias, and is able to evaluate all of them and also can give a summary statistics (MAC Passing Rate) on the sample population.

This work applied and improved the methodology from Rollins et al. ¹⁶. 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 modeling BGC, and provided powerful tools for statistical learning in real life scenario. Nonetheless, these are issues that can be addressed in future research.

Future work will involve running clinical studies under the protocol that subjects will follow when wearing the device such as time stamping for meal size and using only their glucose meter to collect data. If these studies are successful, we plan to develop a prototype armband and evaluate it on several subjects. We envision this device collecting input and output data into the armband where the model will reside. After a sufficient number of lancet measurements have been collected, the model with be built from these data automatically for calibration of the device. After successful calibration, the armband will collect input data, infrequent output data, and display BGC continuously over time on a watch type display or smart phone. Transmission of data from the armband to the display monitor may utilize Bluetooth technology.

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 multiple inputs, subject-specific, modeling under infrequent sampling. However, as this work is only preliminary, there are still several challenges to 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.