|Home | About | Journals | Submit | Contact Us | Français|
Continuous glucose monitors (CGMs) collect a detailed time series of consecutive observations of the underlying process of glucose fluctuations. To some extent, however, the high temporal resolution of the data is accompanied by increased probability of error in any single data point. Due to both physiological and technical reasons, the structure of these errors is complex and their analysis is not straightforward. In this article, we describe some of the methods needed to obtain a description of the sensor error that is detailed enough for simulation.
Data were provided by Abbott Diabetes Care and included two data sets collected by the FreeStyle Navigator™ CGM: The first set consisted of 1032 time series of glucose readings from 136 patients with type 1 diabetes and parallel time series of reference blood glucose (BG) collected via self-monitoring at irregular intervals. The average duration of a time series was 5 days; the total number of sensor-reference data pairs was approximately 20,600. The second data set consisted of 56 time series of glucose readings from 28 patients with type 1 diabetes and a parallel time series of reference BG measured via the YSI 2300 Stat Plus™ analyzer every 15 minutes. The average duration of a time series was 2 days; the total number of sensor-reference data pairs was approximately 7000.
Three sets of results are discussed: analysis of sensor errors with respect to the BG rate of change, mathematical modeling of sensor error patterns and distribution, and computer simulation of sensor errors:
CGM accuracy was modeled via diffusion and additive ARMA noise, which allowed for designing a computer simulator of sensor errors. The simulator, a component of a larger simulation platform approved by the Food and Drug Administration in January 2008 for pre-clinical testing of closed-loop strategies, has been successfully applied to in silico testing of closed-loop control algorithms, resulting in an investigational device exemption for closed-loop trials at the University of Virginia.
Continuous glucose monitors (CGMs) provide detailed time series of consecutive observations on the underlying process of glucose fluctuations. The feedback of such detailed information to patients with diabetes has been shown to have positive influence on their glycemic control, including reduction in glucose variability, time spent in nocturnal hypoglycemia, time spent in hyperglycemia, and levels of glycosylated hemoglobin.1–4 However, a number of studies have concluded that, despite eight years of development, the CGM technology continues to face challenges in terms of sensitivity, stability, calibration, and the physiological time lag between blood glucose (BG) and interstitial glucose (IG) concentration.5–11 Thus, it is frequently concluded that the abundance of information about glucose fluctuations carried by the CGM data stream is to some extent offset by the possibility of sensor errors that exceed in magnitude the errors of the traditional self-monitored blood glucose (SMBG) devices. Such a conclusion, however, is only partially accurate: while the observed error in an isolated CGM data point is indeed generally larger than the error observed in a SMBG data point, the additional information provided by CGM time series allows the application of error-reduction techniques that are unavailable in SMBG devices. For example, deconvolution and other modeling techniques allow for mitigation of certain sensor deviations due to blood-to-interstitial time delay.12,13
The key to CGM error mitigation is detailed analysis and subsequent mathematical modeling, which allow the understanding of the sources and magnitude of sensor errors. The analytical approach proposed in this manuscript is based on two principles: (i) CGMs assess BG fluctuations indirectly—by measuring the concentration of IG—but are calibrated via self-monitoring to approximate BG; and (ii) CGM data reflect an underlying process in time and therefore are a time series consisting of ordered, in-time, highly interdependent data points. The first principle stipulates that calibration errors would be responsible for a portion of the sensor deviation from reference BG.13 Thus, in accuracy studies, the first step of the analysis should be the investigation of calibration errors via simulated recalibration, using all available reference data points. Furthermore, because CGMs operate in the interstitial compartment, which is presumably related to blood via diffusion across the capillary wall, the second step of modeling sensor deviations from BG should be the description of this diffusion process. Models of blood-to-interstitial glucose transport have been proposed and are reasonably well accepted by the scientific community as an approximation of the possible physiological time lag between BG and IG concentration.14–16 The second principle stipulates that the temporal structure of CGM data is important and should be taken into account by the analysis of CGM errors. In particular, established accuracy measures, such as mean absolute/relative difference, present an incomplete picture of sensor accuracy because these measures judge the proximity between sensor and reference BG at isolated points in time, without taking into account the temporal structure of the data. In other words, a random reshuffling of the sensor-reference data pairs in time will not change these accuracy estimates. Thus, in order to account for the dynamics of sensor errors, higher-order temporal properties need to be investigated. A wide array of modeling techniques is offered by time-series methods, such as autoregression, autocorrelation, and spectral analysis. In this paper we propose the use of an autoregressive moving average (ARMA) model to account for the time dependence of consecutive sensor errors.
Finally, the detailed understanding and modeling of sensor errors allows the next step, their computer simulation. This in turn allows the development of a simulated “sensor”, which is useful for in silico testing of diabetes treatment strategies, such as open- or closed- loop control, under the realistic conditions of imperfect CGM.
To decompose the sensor errors, we used techniques from linear regression, kernel density estimation, derivative estimation, and time series analysis, each allowing us to access specific characteristics of the sensor/BG discrepancy. We also provided examples of each analysis using data provided by Abbott Diabetes Care (Alameda, CA).
The data used as an example in this paper comes from two different data sets provided by Abbott Diabetes Care:
To our knowledge, no model-free calibration procedures are available. Therefore, we used the method generally used by sensor producers, i.e., linear/quadratic fitting with time-delay compensation. A posteriori recalibration of the data was performed for each sensor. The major difference between a posteriori and real-time calibration is the availability of all reference BG points; therefore, for a fixed calibration function (the relation between sensor current and BG), a posteriori calibration is optimal in minimizing the sensor readings-reference glucose discrepancy. In this study, a posteriori calibration was performed by linear regression, matching the interpolated sensor readings (if the timing of the reference fell between two readings of the sensor) to the reference measurements. The sum of squares was assessed only at the points of reference measurement. The result of the linear regression was considered to be the recalibrated sensor trace.
Once the sensors were recalibrated, we used kernel density estimation to approximate the distribution of the sensor readings for different glucose references. Each sensor/reference pair is associated with a Gaussian kernel [see Equation (1)] centered about the pair and of predetermined width. More details on the selection of the width and kernel function can be found in Hastie et al.17 The density was then computed as the weighted sum of all kernels:
Here N is the number of pairs, si is the sensor reading of pair i, ri is the reference measure of pair i, and σ is the kernel width. This estimation of the density is slightly biased by the fact that negative glucose values do not occur; the bias is reduced by choosing σ to be less than 25% of the smallest reference value.
Therefore the mean sensor reading (s) and the mean error (ε) for each reference BG was computed using Equation (2).
The premise behind our methodology for delay estimation is that if the CGM measurements are delayed, then the sensor error will be dependent on the rate of change of glucose. Therefore, by computing D(s,r) for different ranges of rates of change in glucose (bins), and studying the differences between bins, we can study the delay. To do so we needed to estimate the density in particular bins using Equation (3):
To compute a robust rate of change we used sliding linear regression computed for consecutive 20-minute windows.
As discussed above, CGMs measure glucose in the interstitial fluid while being used to assess BG, which creates a delay. Therefore, the distribution of sensor errors could not be computed directly by using the difference between reference and sensors at the same time points. Since the physiological delay described above is not constant across subjects or within a subject over a long period of time, we synchronized sensor and reference BG using a first-order diffusion model in the calibration equation [Equation (4)].
Here GI and GB are the interstitial and the blood glucose concentrations, respectively; τ is the diffusion time constant, and α and β are the calibration parameters. Spanning the possible values of τ, we applied the same linear regression technique as before to estimate α and β. The final solution is the set (τ, α, β) which produces the smallest sum of squares. We do not claim that this procedure, derived from the study by Stein et al.,16 models perfectly the transport of glucose from blood to interstitium nor the functional relationship between electrical current in the sensor and glucose values. However, in the absence of an identifiable model of such a process, we followed the parsimony principle in choosing the simplest available one, i.e., a linear transformation added to a first order, gain 1, diffusion.
Once sensor and reference data were synchronized, we computed the differences between sensor and reference and estimates, and the first four central moments of their empirical distribution: mean, variance, skewness, and kurtosis.
Classical time-series techniques were applied to the recalibrated and synchronized sensor signal to determine the time dependency of sensor errors: the autocorrelation function and the partial autocorrelation functions.
The autocorrelation function (ACF) is fairly straightforward: it is computed as the correlation of an error at time t, with the errors at time t + h, where h = nT, n is an integer, and T is a fixed time interval. (Generally T is set to the time difference between reference measures.) Under weak stationary conditions (the mean and variance of the error do not depend on time) the ACF is only dependent on the lag (h) and not on t, and can be computed using Equation (5).
The partial autocorrelation function (PACF) can be best described as the correlation between errors at time t and t + h, h = nT, excluding information transmitted through t + T, t + 2T, t + 3T,…, t + (n − 1)T. It is similar to the concept of best linear predictor and is commonly computed using the Durbin-Levinson algorithm.20 For more details on PACF please refer to Brockwell et al.22
Using data set 1 we studied the sensor response at different reference glucose levels by estimating the probability distribution of the reference/sensor pairs (recalibrated but not synchronized). The distribution is presented in Figure 1, where blue depicts a very low and red a very high probability of occurrence. We observed that sensors tend to read low at high reference values (the reference/sensor pair tends to fall below the diagonal when the reference is above 200 mg/dl) and high at low reference levels (the reference/sensor pair tends to fall above the diagonal when the reference is below 110 mg/dl). Also, the spread of reference/sensor pairs is positively correlated with the reference level: the distribution is flatter at high glucose levels compared to low glucose levels.
As presented in the introduction, it is widely believed that there is a delay between BG and IG. To verify this claim, we applied the same technique described in the previous section but we clustered the reference/sensor pairs by glucose rate of change: 8 bins of rate of change were selected, the number of pairs in each bin was not constant but was always greater than 1000, the number of pairs per bin was also fairly symmetric around a 0 rate of change, i.e. there were roughly as many pairs between -1.5 and -1 mg/dl/min as between 1 and 1.5 mg/dl/min. This distribution of pairs by rate of change corresponds to the usually accepted distribution of rate of change in the field,23 therefore indicating an absence of bias. The distributions for each are presented in Figure 2.
Observing the distributions in Figure 2, particularly the most likely reference/sensor area (red zone in each distribution), we saw that: (i) at a negative rate of change, the sensors tends to read high (red zone above the diagonal); (ii) at a positive rate, the sensors tend to read low (red zone below the diagonal); and (iii) the extent to which the sensor systematically reads high or low is correlated to the amplitude of the rate of change (e.g., the red zone is further above the 45° line in the top left distribution than the bottom left distribution). To verify the last observation, we computed the average reference/sensor discrepancy as a function of the reference glucose for each glucose rate of change bin (j) using the distribution Dj(s,r). The results of this analysis are presented in Figure 3A. We conclude that there exists a correlation between rate of change and average discrepancy, regardless of reference BG levels. Finally, computing the average reference/sensor discrepancy for a specific rate of change across reference values, we compared these averages with the average rate of change in each bin. The results of this analysis are presented in Figure 3B. The average discrepancy is linearly related to the rate of change (R2 = 0.995), and the slope of this linear relation (without offset) gives an estimate of the delay, which in this data set is 17 minutes.
In this section we study the third component of sensor difference from reference, which remains after recalibration and synchronization of sensor and reference data. We refer to this difference as the sensor error. To study more precisely the sensor error and its time dependency, we used data set 2, which contains YSI glucose measurements taken at 15-minute intervals.
The histogram of the sensor error is presented in Figure 4. The sensor error has a mean of 0.76 mg/dl and a standard deviation of 11 mg/dl, but its skewness and kurtosis show that the distribution of the error is not normal. To estimate its distribution we used the Johnson family and obtain the parameters in Table 1.
To characterize the time dependency of the sensor error we computed the ACF and PACF of the sensor error as described in the Methods section. Since the reference glucose measures were equally spaced at 15 minutes, the lags are integer multiples of 15 minutes. The results of this analysis are presented in Figure 5. Using a significance bound based on the 95% confidence interval for Gaussian white noise processes; we conclude that sensor noise is highly correlated across time up to several hours. The lag of autocorrelation depends on each sensor (Figure 5A), but, in general, sensor noise is best predicted by a linear combination of the sensor noise 15 minutes earlier and a random white noise term (PACF cutting off after lag 1 in Figure 5B).
We found that the discrepancy between sensor and reference glucose differ from random noise by having substantial time-lag dependence and other non– independent identically distributed (iid) characteristics (i.e., the error is independent of previous errors and drawn from the same time-independent probability distribution). The components of the discrepancy are therefore modeled as follows:
These equations were included in a glucose homeostasis simulator,22 implemented in Matlab Simulink®. An example is presented in Figure 6, along with validation of the modeling by comparing the sensor error distributions (empirical and from simulation), as well as the partial autocorrelation functions. No difference between empirical and simulated characteristics was apparent. Moreover, a χ2 test showed that no significant difference exists between the observed and the simulated distribution of sensor errors and the distribution of errors of the FreeStyle Navigator (p > .46).
In this paper, we presented a sequence of analytical techniques beginning with the decomposition of CGM sensor errors into errors due to calibration, blood-tointerstitial glucose transport, and additive noise, leading to computer simulation and building of an in silico sensor based on diffusion and time series modeling. We should emphasize that the term sensor “error”, although well accepted, may not be entirely accurate because part of the “error” is due to physiological processes that result in natural deviations of sensor readings from reference BG levels.
Regardless of terminology, sensor deviations from reference BG levels are generally explained by the specifics of CGM technology—the measurement of glucose takes place in the interstitial compartment, which is associated with, but different from blood. However, CGMs are calibrated against BG, and their accuracy is assessed by comparing IG readings to reference BG. This alone creates two sources of CGM errors, or deviations: calibration and the BG-to-IG gradient.
During analysis, the errors due to calibration can be separated by recalibrating the sensor trace using all available reference BG readings, which is the first step of sensor accuracy analysis. The differences between IG and BG are most evident during rapid BG excursions: sensor errors tend to be negative (high CGM readings) when the BG rate of change is negative and positive (low CGM readings) when the BG rate of change is positive, which is potentially explained by the time lag between the two compartments. This second component of sensor error due to the BG-to-IG gradient can be evaluated by a diffusion model—a commonly accepted technique explaining the possible physiological time lag between blood and IG concentration.14–16 After recalibration and accounting for the BG-to-IG gradient, the residual sensor noise appears to be non-white (non-Gaussian) and the consecutive sensor errors are highly interdependent. Thus, the presented analysis of residual (after calibration and diffusion) sensor errors is based on time-series approaches using ARMA noise to account for the interdependence between consecutive sensor errors. As an example, we presented the autocorrelation and the distribution of the errors of the FreeStyle Navigator.
Based on the presented sequence of mathematical models, the computer simulation of sensor errors becomes possible. Such an in silico sensor includes both generic deviations due to physiology and sensor-specific errors due to particular sensor engineering. The generic deviations can be approximated reasonably well by a diffusion model. The sensor-specific errors appear to be generally limited to the order of the ARMA model and the shape of the sensor error distribution. Thus, both types of errors can be described by a few model parameters and included in a simulation. A specific type of error is omitted by this methodology: high frequency errors (period of 1 to 15 minutes) cannot be modeled with our techniques. At this time no data are available to extract the characteristic of the error at that fine a sample frequency (we would need reference values every minute in a large data set). Nonetheless, the development of a simulated sensor becomes feasible and became a component of the larger metabolic system simulator developed at the University of Virginia and the University of Padova. This simulator was accepted by the Food and Drug Administration for pre-clinical testing of closed-loop glucose control strategies. Such in silico testing is now employed by the Juvenile Diabetes Research Foundation Artificial Pancreas Consortium, provides insights on the effectiveness of control algorithms during realistic conditions, and allows the direct transition from in silico testing to clinical trials.
The authors thank Abbott Diabetes Care (Alameda, California) for sharing their database.
The University of Virginia study was supported by the Juvenile Diabetes Research Foundation Artificial Pancreas Project at the University of Virginia.