Autoregressive Modeling of Drift and Random Error to Characterize a Continuous Intravascular Glucose Monitoring Sensor

Abstract

Background:

Continuous glucose monitoring (CGM) devices have been effective in managing diabetes and offer potential benefits for use in the intensive care unit (ICU). Use of CGM devices in the ICU has been limited, primarily due to the higher point accuracy errors over currently used traditional intermittent blood glucose (BG) measures. General models of CGM errors, including drift and random errors, are lacking, but would enable better design of protocols to utilize these devices. This article presents an autoregressive (AR) based modeling method that separately characterizes the drift and random noise of the GlySure CGM sensor (GlySure Limited, Oxfordshire, UK).

Methods:

Clinical sensor data (n = 33) and reference measurements were used to generate 2 AR models to describe sensor drift and noise. These models were used to generate 100 Monte Carlo simulations based on reference blood glucose measurements. These were then compared to the original CGM clinical data using mean absolute relative difference (MARD) and a Trend Compass.

Results:

The point accuracy MARD was very similar between simulated and clinical data (9.6% vs 9.9%). A Trend Compass was used to assess trend accuracy, and found simulated and clinical sensor profiles were similar (simulated trend index 11.4° vs clinical trend index 10.9°).

Conclusion:

The model and method accurately represents cohort sensor behavior over patients, providing a general modeling approach to any such sensor by separately characterizing each type of error that can arise in the data. Overall, it enables better protocol design based on accurate expected CGM sensor behavior, as well as enabling the analysis of what level of each type of sensor error would be necessary to obtain desired glycemic control safety and performance with a given protocol.

Keywords

blood glucose continuous glucose monitoring modeling sensor simulation

Intensive care unit (ICU) patients commonly experience stress induced insulin resistance resulting in hyperglycemia.^1-7 Hyperglycemia is associated with increased morbidity and mortality.^1,3,6,8 Some studies have shown glycemic control (GC) reduces hyperglycemia and improves outcomes.^2,9,10 However, other studies failed to replicate these results,^11-13 with one reporting an increase in mortality.¹³

A critical confounding factor in evaluating these results is the increased incidence of hypoglycemia observed across almost all studies.¹⁴ Hypoglycemia has been independently associated with increased mortality,^15-17 with one study showing increased mortality rates after a single moderate hypoglycemic event.¹⁶ Thus, it is critical that GC protocols treat hyperglycemia safely, as well as effectively, to achieve the potential outcome benefits.¹⁸

Hypoglycemia in GC often occurs where infrequent blood glucose (BG) measurements combine with, often rapid, changes in patient condition and response to care.^19-22 Continuous glucose monitoring (CGM) devices could help provide safety from hypoglycemia, reduce workload by increasing automation, and thus improve GC by providing real-time BG levels. In particular, the increased temporal resolution CGM devices provide can monitor real-time BG trends, allowing more rapid treatment response to highly dynamic changes in patient condition²³ to modify insulin delivery and avoid hypoglycemia.

CGM devices can also reduce the well reported nursing related GC workload, providing more bedside data with lower blood sampling requirements^22,24,25 and improved ergonomics in GC.¹⁹ However, the increased temporal measurement resolution CGM devices provide is still somewhat outweighed by the larger point accuracy errors in these devices due to sensor drift, bias, and random noise.^22,26-31 Trend accuracy is also an important factor, particularly where alarms indicating hypo and hyperglycemia are concerned.³² There is thus a need to model and account for sensor error, preferably in a generalizable modeling method, which would in turn enable optimal (model-based) design of GC protocols maximizing CGM advantages and minimizing their disadvantages, as well as matching recent consensus statements from medical and industry based working groups.³³

Three CGM models have been developed in the past, primarily for interstitial CGM devices, by Breton and Kovatchev,³⁰ Lunn et al,³¹ and Facchinetti et al.²⁶ The model developed by Breton and Kovatchev has been used for in silico preclinical trials to simulate the effectiveness of using an interstitial CGM device for closed-loop control.^34,35 However, Facchinetti et al have since shown that the modeling methodology used by Breton and Kovatchev may be sensitive to small errors in CGM data calibration or errors in the description of BG-to-interstitial glucose (IG) dynamics.³⁶ Lunn et al³¹ developed a more refined version of the model developed by Breton and Kovatchev by fitting a dynamic model with forcing functions. However, it suffers from the same issue as Breton and Kovatchev’s model. In addition, neither model separately considered sensor drift, thus including it in point error, and thus were not as accurate for sensors where drift occurs.

Facchinetti further developed a model of sensor error incorporating BG-to-IG dynamics, using an autoregressive model to account for additive measurement noise and a linear time-varying model to account for calibration and sensor drift.²⁶ Again, BG-to-IG dynamics are added because the CGM device modeled was an interstitial device. This modeling method was able to account for sensor drift, and also split sensor error into multiple components, including error arising from calibration and measurement noise. However, this modeling method is very data- and labor-intensive, requiring multiple CGM devices per patient and 15 minute intermittent BG measurement intervals. This large amount of data may not be available in past data acquired from CGM device trials, while the large workload required to measure BG every 15 minutes may be a barrier for further sensor error characterization, and thus to modeling and simulation of sensor behavior in new CGM devices.

In addition, none of the previous modeling efforts have taken into account trend accuracy of the CGM devices studied. Signal et al³² developed the Trend Compass and Trend Index to assess a CGM sensor’s trend accuracy, which could be seen as equally important as the measurement of mean absolute relative difference (MARD) to assess point accuracy in the previous modeling efforts.

The GlySure (GlySure Limited, Oxfordshire, UK) CGM device considered in this work is from a newer class of ICU devices measuring venous BG via an intravenous line, thus avoiding IG dynamics. The sensor is comprised of microporous and dialysis membrane, hydrogel, optical fiber and a thermocouple, while the glucose detecting chemistry used is a fluorescent diboronic acid receptor, embedded within the hydrogel. Placement of the sensor itself can occur through either a central venous catheter, or a radial artery catheter.³⁷ This article presents a novel autoregressive (AR) method and model characterizing this CGM sensor. The modeling method is capable of characterizing the CGM sensor with less data than required by previous characterization methods in terms of sensors per patient, and explicitly accounts for sensor drift, while maintaining both the point and trend accuracy of simulated sensor errors. This model is developed and compared to clinical data to assess its validity. The overall modeling method and approach is generalizable to similar devices.

Methods

Sensor modeling

Clinical Data

Data were sourced from an observational pilot trial of the CGM device on 33 cardiac intensive care patients (duration 21-51 hours per patient), where CGM readings were not used clinically for GC. The sensor provides a new reading 4 times per minute, and intermittent BG measures are used to calibrate the sensor (calibration or recalibration BG) approximately every 8 hours. Intermittent independent BG measures not used to calibrate the sensor (reference BG) were taken approximately every 2.5 hours. Each intermittent BG measurement was taken using YSI 2300 STAT Plus (Yellow Springs Instruments, Yellow Springs, OH) or the i-STAT (Abbott Laboratories, Abbott Park, IL), which are highly accurate measures, to minimize the error in reference BG values,^38,39 as reported in Crane et al.³⁷

Patient details and details of the pilot trial can be found in Table 1 and Table 2 from.³⁷ Data from patients 10 and 22 were later discarded due to sensor failure. MARD was calculated between paired BG (calibration and reference) and CGM measurements. The average global MARD (excluding patients 10 and 22) was 9.6%.

Table 1.

Patient Details.

	Cardiac patients
n	33
Duration (hours) (mean) (range)	40.8 (21.1-50.7)
Male (n) (%)	22 (66.7%)
Female (n) (%)	11 (33.3%)
Individuals with diagnosed diabetes (n) (%)	14 (42.4%)
Hypertensive (n) (%)	15 (45.5%)
BMI (mean) (range)	25.3 (17.7-35.8)
Age (mean) (range)	50.8 (19-77)

Table 2.

Pilot Trial Details.

Patient number	No. hours of CGM	No. calibration BG	No. reference BG	MARD	Median [IQR] sensor error (%)
1	48.33	3	14	8.62	8.4 [4.6, 12.3]
2	40.97	3	13	6.20	5.1 [1.7, 10.4]
3	42.66	3	13	6.00	2.9 [1.7, 9.4]
4	47.33	2	15	10.66	9.7 [6.3, 14.3]
5	47.05	3	15	6.85	4.6 [2.0, 11.2]
6	50.69	3	14	9.57	9.3 [7.8, 10.8]
7	46.90	3	14	12.56	13.3 [6.5, 16.8]
8	43.43	3	14	4.47	2.7 [1.7, 5.9]
9	47.05	7	14	13.25	10.9 [3.5, 18.2]
10	—	—	—	—	—
11	48.69	3	15	15.62	8.5 [4.1, 11.2]
12	42.84	3	13	11.73	6.1 [1.5, 12.0]
13	40.91	3	12	6.10	14.5 [10.8, 25.4]
14	45.60	4	13	12.10	13.0 [8.1, 21.3]
15	39.37	5	10	4.60	12.8 [6.0, 18.3]
16	43.51	7	9	6.43	6.9 [2.7, 13.7]
17	37.71	3	12	7.48	10.3 [6.1, 21.5]
18	36.72	5	11	7.84	11.3 [9.1, 15.7]
19	39.23	2	13	4.85	4.5 [2.7, 8.5]
20	40.42	4	11	7.19	4.0 [2.3, 7.2]
21	36.35	4	11	8.27	7.2 [4.0, 7.9]
22	—	—	—	—	—
23	37.77	3	13	5.29	6.0 [2.7, 10.3]
24	39.03	4	12	12.20	3.5 [1.5, 8.0]
25	38.12	3	13	4.77	12.8 [5.0, 17.4]
26	36.69	3	12	6.22	4.4 [2.2, 6.4]
27	40.12	4	12	12.36	4.7 [3.2, 7.7]
28	37.65	4	12	14.49	8.5 [7.7, 16.0]
29	37.81	3	12	22.93	16.5 [9.7, 20.9]
30	36.51	4	11	10.12	21.9 [18.5, 27.7]
31	37.54	4	12	8.19	8.3 [3.3, 13.9]
32	21.14	2	7	17.38	7.0 [3.4, 11.8]
33	25.49	3	9	14.35	14.0 [11.9, 17.7]

Patients 10 and 22 were discarded in this retrospective analysis due to sensor failure.

Sensor Characterization

Sensor characterization uses two independently defined AR models to separately capture drift and higher frequency sensor fluctuations, where most other methods have not explicitly accounted for sensor drift.^30,31 This method is able to be used on the clinical data where there was only one CGM sensor per patient, which would not have been possible with the method of,²⁶ as the characterization of separate noise components requires multiple sensors for each patient.

Drift

Clinical data are divided into separate periods between recalibration points for each patient. Drift is characterized for any given patient trace between recalibration measurements using the percentage difference between sensor and reference measurements, as assessed half-hourly using sensor glucose (SG) ( $B G_{S G / 30 m i n}$ ), and intermittent BG interpolated ( $B G_{I M / 30 m i n})$ between calibration and reference measurements ( $B G_{I M})$ . Interpolated intermittent BG_IM measures are a vector with an interpolated value every 30 minutes, defined:

B G_{I M / 30 m i n} = {i n t e r p (B G_{I M}) |}_{t = 0 : 30 : t_{e n d}}

Half hourly sampling of the CGM sensor trace simply takes the paired CGM value at that time:

B G_{S G / 30 m i n} = {s a m p l e (B G_{s e n s o r}) |}_{t = 0 : 30 : t_{e n d}}

The half hourly interval matches observed physiological and clinical time frames for BG fluctuations and trends, and thus captures the broad trends caused by sensor drift and, over several samples, eliminates the impact of random errors at any given point. Thus, a potential percentage drift can be calculated at each paired half hourly value:

D r i f t = \frac{B G_{S G / 30 m i n} - B G_{I M / 30 m i n}}{B G_{I M / 30 m i n}}

A lag-2 AR model is then used to characterize the observed drift for a given patient’s data or for the whole cohort. This AR model uses the entire cohort’s data (N = 31 patients), and is defined:

D r i f t_{n + 1} = α_{d} + β_{d} * D r i f t_{n - 1} + γ_{d} * D r i f t_{n} + ξ_{d}

Model parameters α_d, β_d and γ_d are identified using linear least squares from the half hourly [Drift_n+1, Drift_n, and Drift_n-1] data points derived from the entire clinical data cohort (N = 31 patients) and Equations 1 -3, assuming ξ_d = 0. The parameters α_d, β_d and γ_d thus capture sensor drift behavior over the entire cohort.

Having identified the best fit (α_d, β_d and γ_d) for the cohort, a drift noise term is calculated by rearranging Equation 4 and solving for the residuals, ξ_d. Outlier drifts from Equation 3 of more than ±25%, where 99.3% of the data are within this threshold, were discarded to maintain sufficient data density for probability modeling. These results are used to create a drift noise model by smoothing a continuous distribution function across the ξ_d residual results obtained by rearranging Equation 4.

Random Sensor Noise

Random sensor noise fluctuations are assessed at the CGM sensor sample rate (4 measures per minute) and are sampled from the interpolated BG and CGM device measurements every 15 seconds or 0.25 minutes, yielding:

B G_{b a s e / 0.25 m i n u t e s} = {i n t e r p (B G_{S G / 30 m i n}) |}_{t = 0 : 0.25 : t_{e n d}}

The fractional difference between these linearly interpolated ‘base’ BG points and the real sensor trace is defined:

S e n s o r F l u x = \frac{B G_{s e n s o r} - B G_{b a s e / 0.25 m i n u t e s}}{B G_{b a s e / 0.25 m i n u t e s}}

Another lag-2 AR model is used to characterize these sensor fluctuations around the trend, defined:

\begin{array}{l} S e n s o r F l u x_{n + 1} = α_{s f} + β_{s f} * S e n s o r F l u x_{n - 1} \\ + γ_{s f} * S e n s o r F l u x_{n} + ξ_{s f} \end{array}

The model parameters α_sf, β_sf and γ_sf are identified using the entire data cohort (N = 31 patients) and linear least squares from [SensorFlux_n+1, SensorFlux_n, and SensorFlux_n-1] data points obtained every 15 seconds (0.25 minutes), as derived from clinical data, and Equation 6, assuming ξ_sf = 0.

After identifying a linear best fit (α_sf, β_sf and γ_sf) across the entire data cohort, the sensor fluctuation noise term is calculated by rearranging Equation 7 and solving for the residuals, ξ_sf. Outlier fluctuations of more than ±1%, where 99.9% of the data are within this threshold, were discarded. These results are used to create a sensor fluctuation noise model, similar to the drift noise model, by smoothing a continuous distribution function across the data range of ξ_sf.

Illustrated Example of Sensor Characterization

An example of the modeling process is shown using data from patient 2 in Figure 1. Figure 1a shows the clinical sensor data for patient 2 over the first 16 hours, with the data split into separate periods between recalibration times to characterize drift. Figures 1b and 1c show characterization of the AR drift model using Equations 3 and 4, where Figure 1c shows the plane of best fit for model parameters α_d, β_d and γ_d, and the drift noise term definition (ξ_d). Figures 1d-1fshow characterization of the sensor fluctuation model using Equations 6 and 7, with planes of best fit for model parameters α_sf, β_sf and γ_sf and the sensor fluctuations noise term definition (ξ_sf).

Figure 1.

Steps of the sensor characterization process, showing how the drift and sensor fluctuations were characterized. (a) Clinical data plotted for patient 2. (b) Patient 2 with BG resampled every half hour. Red arrows show how drift was characterized (Drift = (BG_SG/30min-BG_IM/30min)/(BG_IM/30min). (c) Drift_n+1 plotted against current (Drift_n) and previous (Drift_n-1) for all patients. The plane of best fit is shown in black, with a red arrow showing how the residual ξ_d is sampled. (d) Patient 2 with BG resampled every minute. Red arrows show how the sensor fluctuation was characterized (Note small x and y scales). (e) Sensor fluctuation_n+1 plotted against current sensor fluctuation (n) and previous sensor fluctuation (n-1) for all patients. The plane of best fit is shown. (f) Sensor fluctuation_n+1 plotted against current sensor fluctuation (n) and previous sensor fluctuation (n-1), with a red arrow showing how the ξ_sf is sampled for the fluctuations.

Sensor Simulation

Sensor simulation for validation simulates a sensor trace given intermittent BG_IM measures in a process essentially the reverse of sensor characterization. Intermittent BG_IM measures are used as a base to simulate the CGM sensor, as that is all the clinical data that might be available when simulating a virtual patient to be monitored by CGM.^40-43 They are interpolated half hourly, and drift applied using Equation 4. The resulting BG with drift is defined:

B G_{b a s e (t = 0 : 30 : t_{e n d})} = B G_{I M / 30 m i n u t e s} * (1 + D r i f t)

Where Drift is calculated using Equation 4, with ξ_d drawn randomly from the cohort probability distribution generated from the clinical data. At t = 0 and any calibration BG, approximately every 8 hours for this device, the condition [Drift_n, Drift_n-1] = [0.0, 0.0] is used to recalibrate the CGM sensor to match the simulated intermittent BG measurement, providing a point to point calibration. This approach could be modified to account for any more complex calibration process.

BG_base is then linearly interpolated to provide enough data points to match this CGM device’s 4 times per minute rate, such that:

B G_{b a s e / 0.25 m i n u t e s} = {i n t e r p (B G_{b a s e}) |}_{t = 0 : 0.25 : t_{e n d}}

The simulated sensor output is then defined:

\begin{array}{l} B G_{n e w_s e n s o r (t = 0 : 0.25 : t_{e n d})} = B G_{b a s e / 0.25 m i n u t e s} \\ * (1 + S e n s o r F l u x) \end{array}

Where SensorFlux is calculated according to Equation 7, with ξ_sf drawn randomly from the cohort probability distribution generated from clinical data.

Once again, at t = 0 and any calibration BG, the condition $[S e n s o r F l u x_{n}, S e n s o r F l u x_{n - 1}] = [0.0, 0.0]$ is used to recalibrate the CGM sensor. Finally, an additional limit of a maximum drift of 40% was applied to the drift AR model matching extremes seen in the clinical sensor data.

Regarding calibration, setting $[D r i f t_{n}, D r i f t_{n - 1}] = [0.0, 0.0]$ and $[S e n s o r F l u x_{n}, S e n s o r F l u x_{n - 1}] = [0.0, 0.0]$ effectively applies a point to point recalibration, allowing the sensor trace to go through a recalibration point. Divergence from the interpolated BG is then initiated by the ξ_d and ξ_sf terms. Figure 2 outlines the steps during simulation, with a sensor trace and virtual patient forward simulated for the first 4 hours given the initial intermittent BG data in Figure 2a.

Figure 2.

Sensor simulation steps, calibration not shown. (a) An example patient. Intermittent clinical BG measures are interpolated. (b) AR drift is applied to the interpolated BG at half hour intervals (Equation 8). (c) AR sensor fluctuations are added to the interpolated AR drift (Equations 9 and 10). (d) Comparison of intermittent BG and final sensor simulation.

Sensor Model Validation: Qualitative and Quantitative

Sensor traces are simulated using intermittent BG from the clinical data cohort from which the model was built. To test consistency between the model and the clinical sensor data, several sensor simulations were overlaid with the clinical data for each patient trace and compared in a blinded test. If the clinical sensor traces were difficult to visually distinguish from simulation, then the model was qualitatively accepted as broadly capturing key behavior. This qualitative assessment enables assessment of trends and features not easily compared in quantitative tests.

To quantitatively validate the model, a single simulation was run to generate a single sensor trace from the sensor model for each patient to compare to the clinical sensor data. Percentage difference distributions for the simulated sensor data and the clinical sensor data were compared, and a Clarke error grid (CEG) plot and Bland-Altman plot constructed. The Bland-Altman plot enables analysis of any differences in bias behavior between clinically measured and modeled sensor traces. Trend accuracy was assessed and compared for the clinical data and simulation data using the Trend Compass described by Signal et al,³² where a Trend Compass and Trend Index (defined in Equation 11), an absolute mean angle from perfect trend accuracy, were produced for both clinical and simulation data to validate the reproduction of trend accuracy using the simulation method. The Trend Compass is a visual tool that can be used to quickly identify the trend accuracy of a particular CGM sensor. The compass itself is split into four separate quadrants depending on the relative changes of BG and SG, with measurements plotted according to the angle between the interpolated BG vector and the corresponding CGM vector (θ), and radially according to the last interpolated BG measurement. The angle of theta (θ) is plotted as measured from the vertical and represents the level of agreement between the rates-of-change in BG and SG, where points plotted on the vertical have perfect trend accuracy. It thus evaluates rising (upper half) and falling (lower half) trends separately, where the two sides distinguish trends where CGM measured SG rises (or falls) faster or slower than BG. This quantitative validation was carried out for each patient and sensor trace in the clinical cohort.

T I = \frac{1}{n} \sum_{i = 1}^{n} | θ (i) |

Further Monte Carlo simulations were undertaken, until a total of 100 simulations worth of data were generated for each patient. Individual patient simulation traces were plotted on top of the clinical data to further check the consistency between the clinical data and simulated traces. The range of simulated sensor behavior was also compared to the original clinical data to ensure simulated sensor behavior reflected extremes in actual sensor behavior. For each patient, the minimum simulated BG value was taken for each time point over the 100 simulations to generate an overall minimum sensor profile. A maximum sensor profile was also produced for each patient in a similar method, using the maximum simulated BG values for each time point. Together these minimum and maximum profiles formed an area profile of all possible simulated BG values, which were then compared to the clinical sensor traces for each patient.

Convergence Analysis

To test that this model methodology is able to work with limited reference readings, the number of reference BG measurements that were not recalibration measures was halved (42% less measures), and the sensor was recharacterized for drift and noise. The global MARD for the simulations of the sensor was reevaluated and percentage difference distribution and CEG plots generated to compare with the previous analyses.

Results

Sensor characterization

Table 3 gives sensor model parameters identified from the cohort data, and Figure 3 shows the noise term model distributions (ξ_d and ξ_sf), raw and smoothed fit, for drift and sensor fluctuations.

Table 3.

Autoregressive (AR) Model Parameters for Drift and Sensor Fluctuations.

AR pass	Key characteristic	$α$	$β$	$γ$	$ξ$ median	$ξ$ max (absolute)	R ²
1	Drift	–0.00147	–0.07152	0.9698	–0.0023	0.2486	0.81
2	Sensor fluctuations	–6.0e–6	–0.261	1.261	–2.66e–06	0.01	0.99

Figure 3.

Random noise distributions for autoregressive drift and sensor fluctuations (fraction). (a) Noise model for drift, ξ_d. (b) Noise model for sensor fluctuations, ξ_sf.

Sensor Simulation

Figure 4 shows four example patients with the real sensor trace plotted alongside two simulated traces. Sensor behavior is visually consistent between simulated and real sensor traces, with drift and sensor fluctuations of approximately the same magnitude across all traces. In a few sensor traces, such as the clinical data in Figures 5b, there is additional high frequency low amplitude noise present the model cannot capture. This noise is not likely to affect glycemic control applications of this sensor, whether in simulation or practice. The model is thus considered qualitatively good.

Figure 4.

Clinical sensor data with two sensor traces generated from the model. (a) Clinical sensor data plotted alongside 2 simulations for patient 1. (b) Clinical sensor data plotted alongside 2 simulations for patient 2. (c) Clinical sensor data plotted alongside 2 simulations for patient 4. (d) Clinical sensor data plotted alongside 2 simulations for patient 5.

Figure 5.

Percentage difference distribution plot of the clinical sensor data vs the simulated data.

Figure 5 shows the distribution of percentage differences over all N = 31 patients clinical and simulation data, with one simulation per patient. The percentage differences were calculated by taking SG, subtracting the interpolated intermittent BG, and then dividing through by the interpolated intermittent BG. The model simulated data distribution is slightly tighter than the clinical data. This outcome is mainly due to the point-to-point recalibrations used in the simulations, which is a slightly more accurate recalibration than the one used in the clinical trial. There is also some slight non-Gaussian distribution of percentage differences at the extremes, compared to the Gaussian noise distribution used in the AR modeling. However, the distribution is still very similar, indicating that the modeling method is able to accurately recreate the percentage differences in the measurements of intermittent BG and SG.

The CEG in Figure 6 shows the model behaves in a consistent manner to the clinical data. While the distributions are consistent across the BG range, as expected, there are a few outliers in the clinical data not captured by the model. The percentages of measurements falling within the zones of the CEG plot are shown in Table 4. Slightly more data points fell within zone A than zone B when comparing the simulated results to the clinical results, suggesting that the simulation method may be slightly more accurate than the clinical sensor. However, this difference is minor (~6.6% change between zones A and B), and thus the model could still be considered to quantitatively represent the sensor behavior well. In Figure 7, the Bland-Altman plot shows no significant bias across the observed BG range, and the plotted lines of ±2σ for the clinical and simulated data show the strong similarity between the model outputs and the clinical sensor data over this 95% range.

Figure 6.

Clarke error grid plot of the clinical sensor data and the clinical simulated data, both resampled half hourly.

Table 4.

Percentages of Simulated and Clinical Measurements Falling Within the Zones of the CEG.

Zone	Simulated data (%)	Clinical data (%)
A	89.3	82.7
B	10.6	17.1
C	0	0
D	0.1	0.2
E	0	0

Figure 7.

Bland-Altman plot of clinical sensor data and simulated clinical data.

Of note, there is very little clinical data below 90 mg/dl. Underlying sensor model assumptions apply constant sensor behavior across the full BG range resulting in similar proportional BG error at high and low BG. This assumption is used for lack of other data from the sensor at this time. In this case, at lower BG, this choice translates to a consistent percentage error, resulting in slightly lower absolute BG errors.

Figure 8 shows the Trend Compass³² plot for the clinical sensor data and simulated data. The simulated data matches the clinical data well, indicating that the model is able to capture the trend accuracy of the clinical sensor as well as the point accuracy. The Trend Index, as described in Signal et al,³² of the clinical sensor and simulated sensor were 10.9° and 11.4° respectively, while the IQR of the theta values used to evaluate the Trend Index were [3.4°, 16.2°] and [4.0°, 16.1°] for the clinical data and simulated data, respectively. The similarity between the Trend Indexes and the IQRs of theta values further show that the behavior of the model is consistent with the clinical data, particularly important for trend simulation which has not been covered in other models. Overall, the model quantitatively represents the sensor behavior well.

Figure 8.

Trend Compass plot of clinical sensor data and simulated clinical data.

Table 5 compares the median clinical SG for each patient, as measured by the CGM sensor at each intermittent BG measurement time, to the average median value from the 100 model simulations of SG. The measured median values for each patient are comparable and differ only slightly, as shown by the percentage error, which had median and IQR range values of 1.2%, –1.1%, and 3.1%, respectively. Differences are primarily due to the application of a cohort-model to individual patients, confirming that the model is able to capture the average sensor behavior over the cohort well. The largest percentage error occurred for patient 32, which had a clinical sensor reading that was consistently lower than the intermittent BG measurements. This different CGM sensor behavior could be due to sensor malfunction in the clinical trial, and is not necessarily a failure of the characterization method or simulation method.

Table 5.

Median SG for the Clinical Data and the Mean of the Median Simulated SG Value Over 100 Simulations for Each Patient.

Patient	Median clinical SG (mg/dl)	Average median simulated SG (mg/dl)	Percentage error (%)
1	174.6	167.4	4.1
2	124.2	120.6	2.9
3	169.2	158.4	6.4
4	181.8	180	1.0
5	149.4	147.6	1.2
6	154.8	156.6	–1.2
7	145.8	147.6	–1.2
8	154.8	158.4	–2.3
9	142.2	140.4	1.3
10	—	—	—
11	144	138.6	3.8
12	174.6	172.8	1.0
13	181.8	189	–4.0
14	142.2	147.6	–3.8
15	178.2	181.8	–2.0
16	172.8	172.8	0.0
17	145.8	140.4	3.7
18	174.6	171	2.1
19	174.6	169.2	3.1
20	199.8	205.2	–2.7
21	169.2	163.8	3.2
22	—	—	—
23	183.6	183.6	0.0
24	232.2	225	3.1
25	145.8	140.4	3.7
26	162	163.8	–1.1
27	169.2	167.4	1.1
28	194.4	187.2	3.7
29	151.2	147.6	2.4
30	160.2	154.8	3.4
31	151.2	147.6	2.4
32	171	208.8	–22.1
33	154.8	153	1.2

Patients 10 and 22 were discarded in this retrospective analysis due to sensor failure.

Figure 10 shows examples of the 100% range of the simulated sensor traces plotted against the clinical CGM sensor data for the same 4 patients in Figure 9. These plots show the range of simulated SG that could be simulated for the given “true” actual intermittent BG trajectory experimentally measured and its comparison to what the The average global MARD for each simulated patient trace over the 100 sensor simulations, excluding patients 10 and 22, was 9.6%, where the range of the global MARD of each simulation was from 8.3% to 10.9% over all 31 patients. These values compare well with the clinical global MARD of 9.9%.

Figure 9.

20 sensor simulations (out of 100 total) plotted for some example patients, with clinical reference and calibration BG plotted for reference. (a) Patient 1. (b) Patient 5. (c) Patient 6. (d) Patient 8.

Figure 10.

Clinical CGM sensor data (blue) plotted against the range profile generated with 100 Monte Carlo simulations (gold). (a) Patient 1. (b) Patient 5. (c) Patient 6. (d) Patient 8.

Figure 9 shows examples of four typical patients individual simulations, with first 20 of the 100 total simulations plotted for each patient. They each show the multitude of potential sensor traces for a patient given their particular intermittent BG measurements, with a single (arrow) recalibration point included for each patient. Note that each sensor simulation, regardless of the amount of sensor drift at the particular point in time, collapses to the recalibration BG at the recalibration time from the clinical trial, before initiating divergence again through the ξ terms in the model.

sensor recorded. Importantly, this range is calculated over 100 simulations, but no trace follows, for example, the outlying line as the range shows the width of possibilities at any given point. The clinical data fell within the simulated sensor range for a large majority of time (> 95% for 88% of patients) or all the time (70% of patients) for most patient profiles.

Convergence Analysis

Figure 11 shows the percentage difference distribution plot of the simulated sensor using the recharacterized model of the reduced reference measurement data set. The plot is very similar to the plot of the full sensor model in Figure 5, where the minor differences are able to be explained by the random nature of the CGM sensor trace simulation method. Figure 12 also shows the CEG plot of the recharacterized model simulation with reduced reference measurements. The plot is very similar to the plot of the full sensor model in Figure 6, while Table 6 gives a breakdown of the percentages of measurement pairs falling within the zones of the CEG plot. Again, the percentages are very similar to the percentages of the full sensor model in Table 4, with minor differences also being able to be explained by the random nature of the simulation method. The global MARD for the reduced reference measurement simulation was calculated to be 10.1%, which again is similar to the full sensor model simulation global MARD of 9.6%.

Figure 11.

Percentage difference distribution plot of the recharacterized model simulation with reduced reference measurements.

Figure 12.

CEG plot of the recharacterized model simulation with reduced reference measurements.

Table 6.

Percentages of Simulated and Clinical Measurements Falling Within the Zones of the CEG.

Zone	Simulated data (%)	Clinical data (%)
A	86.1	82.7
B	13.7	17.1
C	0	0
D	0.2	0.2
E	0	0

Discussion

Sensor Model

Overall, a sensor model was developed from clinical data and simulated measured sensor behavior well. The MARD, CEG plot and Trend Compass plot for the simulated sensor and the clinical data are very similar (MARD 9.6% vs 9.9%, Trend Index 11.4° vs 10.9° respectively). Equally, simulated sensor traces were difficult to visually distinguish from clinical data, and clinical data fell within the simulated sensor range for a large majority of time (> 95% for 88% of patients) or all the time (70% of patients) for all patient profiles. Overall, results suggest the method and model are qualitatively and quantitatively able to describe virtually all the observed CGM sensor behavior, including accurately capturing trend behaviors, which is important for future work testing the sensor in simulation for CGM-based protocol development and optimization.

This method of sensor error characterization compares well with the method developed by Facchinetti et al,²⁶ with both methods able to achieve a similarity between the simulated global MARD and the global MARD from clinical studies of the respective devices studied.^26,44,45 The method presented here and the method developed by Facchinetti et al also have the advantage of explicitly and independently characterizing sensor drift compared to other prior models.^22,30,31 The method presented here could also be perceived as more straight-forward to implement than that of the method developed by Facchinetti, requiring much less data, and also could be used on existing data sets where only data from one CGM device per patient is available, which is currently typical.

Another possible advantage of the method presented here is that it may require less reference BG to characterize the error from a sensor and simulate sensor error compared to.²⁶ The clinical data had reference BG measurements every 1-4 hours, which were then interpolated to carry out the characterization step. The clinical data Facchinetti et al used required a reference BG measurement every 15 minutes, possibly to enable a more accurate parameter identification of the parameters used in the model. While this high level of clinical data allows precise sensor characterization, it is highly intensive and costly to gather and would not be available in typical clinical or regulatory trial data. Thus, the model created can be more easily developed and implemented.

The reproduction of trend accuracy within CGM sensor simulation is an important advantage of this modeling method that has yet to be tested on other modeling efforts in other publications. Poorer trend accuracy in simulation than seen in the clinical sensor can result in more episodes of undetected hypoglycemia, or a higher rate of false alarms than what would be seen in practice.³² Conversely, higher trend accuracy in the simulated sensor than in the clinical device would hide the number of missed hypoglycemic events and give a lower rate of false alarms than what would occur in practice, making the simulation performance seem much better than clinical performance. Both scenarios are mitigated by the accurate trend accuracy reproduction by the model presented in this article.

A key result of having a model and generalizable sensor modeling method for a CGM sensor that can accurately reproduce point accuracy and trend accuracy, is being able to test the effects of using CGM sensor readings in place of intermittent BG readings during GC in a virtual environment. The model can then be used in virtual patient trials to optimize a CGM-enabled GC protocol for any characterized CGM sensor. More importantly, such models, and a generalizable method for making them, also enable the ability to assess what level of CGM performance makes the technology feasible in the ICU to safely improve care and reduce workload.

New developments in CGM devices also have resulted in BG being measured at a higher frequency, such as the device used in this study that is able to measure BG 4 times per minute. A benefit of the model presented here is the ability to capture both the fast and slow random error dynamics of the CGM device. The capturing frequency can be adjusted to match the measurement frequency of other CGM devices that are to have the sensor errors characterized and simulated.

Limitations

One of the limitations of this study is the limited BG < 90 mg/dl in the clinical data used to construct the model. As a result, use of this model outside the BG range used to generate it invokes the assumption that sensor behavior is consistent as a percentage across the BG range, further implying smaller absolute errors at lower BG, and higher absolute errors at higher BG. This assumption of sensor behavior is consistent with some studies carried out utilizing brands of interstitial CGM devices,^44,46,47 while others report an increase in MARD at lower BG.^48,49 Further research would need to target data in this BG range to develop more a reliable sensor model at lower BG.

The linear assumption in the interpolation of BG can also be questioned, as glucose levels could rise or fall faster or slower than the expected linear line. This change depends on changes in patient condition, as well as any changes in insulin and nutrition inputs, which are not known to us in this dataset. In such conditions, recent work has shown, including times when the input data are better known, that a linear assumption is the best in terms of overall error distributions.⁵⁰

Another limitation is sensor signal delays. Both the signal filter and the diffusion through the sensor membrane will introduce some delay, although it is likely a smaller signal delay of 10-30 seconds. However, with the relatively high frequency of sensor measurements the sensor makes (4 times per minute) and the relatively long period of times between treatment decisions for which these measurements might be used for in a glycemic control protocol (0.5-3 hours), this combined signal delay would not be a hugely significant source of error clinically, and thus was not explicitly accounted for in the development of the model. Future work could include an analysis of signal delay to fully evaluate this hypothesis and evaluate its effect on performance in glycemic control.

A further limitation arises from the clinical sensor data. There were only 1312 hours of recorded data and the median recording period was 39.4 hours. The data that this model was produced from was also the data that this model was tested on, limiting the conclusions of this analyses. More data, particularly from patients that stay for longer than 37.6 hours would allow accurate modeling of long term sensor behavior, particularly if errors change over time in situ, as occurs with interstitial sensors.^48,51-53 Testing the modeled sensor on other sets of data would also help further reinforce this work.

However, it should also be noted that the convergence analysis testing lesser numbers of measurements and model accuracy showed that reducing the reference measurements by 42% before creating the model did not have not significant effect. In particular, the MARD and CEG results were still very close to those of the device in Figures 5 and 6. Thus, we can also conclude that unless specific dynamics are missing from these patients in this dataset, the number of hours used to create this model is acceptable in its ability to capture the fundamental dynamics

Finally, the analysis was limited by the amount of patient data from the clinical trial. A larger clinical cohort would have allowed more in-depth analysis through cross-validation. In addition, it would enable more precise characterization of drift ranges and inclusion of larger drifts in the model would not have skewed results. However, the results from the MARD, CEG plot, Bland-Altman plot and Trend Compass plot indicate that there was enough data to generate a good model that closely captures over 99% of the observed data, since exclusions in sensor modeling eliminated less than 1% of the clinical data.

Conclusions

A CGM sensor was characterized from patient clinical data using an AR modeling approach. The method presented here has the benefits of explicitly accounting for sensor drift and requiring far fewer independently sampled blood glucose measures than other methods. Sensor traces can be simulated for BG taken at a clinically realistic rate to create the model. Sensor simulations showed modeled sensor behavior was very similar to the original clinical data, with very high similarity in MARD, and equally similar Bland-Altman and CEG results further validating the model. The novel use of the Trend Compass to validate the trend accuracy reproduction within simulation further showed that the model method is able to accurately capture both point accuracy and trend accuracy. The overall model method is general to any similar sensor and readily extended to interstitial sensors, with or without including interstitial glucose dynamics. It is easily simulated on typical clinical data and thus readily able to be incorporated into proven virtual patients to optimize protocol designs to utilize CGMs in the intensive care unit for glycemic control.

Footnotes

Abbreviations

AR, autoregressive; BG, blood glucose; CEG, Clarke error grid; CGM, continuous glucose monitoring; d, drift; GC, glucose control; ICU, intensive care unit; IG, interstitial glucose; IM, intermittent; MARD, mean absolute relative difference; SF, SensorFlux or sensor fluctuations; SG, sensor glucose.

Declaration of Conflicting Interests

The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Funding

The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: Clinical data and sensor information were provided by GlySure Limited (Oxfordshire, UK), though GlySure was not involved in the modeling methodology or patient simulation. The authors also acknowledge the support of the EUFP7 and RSNZ Marie Curie IRSES program, the Health Research Council (HRC) of New Zealand, the MedTech CoRE and TEC, and NZ National Science Challenge 7, Science for Technology and Innovation.

References

Capes

Hunt

Malmberg

Gerstein

HC.

Stress hyperglycaemia and increased risk of death after myocardial infarction in patients with and without diabetes: a systematic overview. Lancet. 2000;355(9206):773-778.

van den Berghe

Wouters

Weekers

et al . Intensive insulin therapy in critically ill patients. N Engl J Med. 2001;345(19):1359-1367.

Umpierrez

Isaacs

Bazargan

You

Thaler

Kitabchi

AE.

Hyperglycemia: an independent marker of in-hospital mortality in patients with undiagnosed diabetes. J Clin Endocrinol Metab. 2002;87(3):978-982.

Mizock

BA.

Alterations in fuel metabolism in critical illness: hyperglycaemia. Best Pract Res Clin Endocrinol Metab. 2001;15(4):533-551.

McCowen

Malhotra

Bistrian

BR.

Stress-induced hyperglycemia. Crit Care Clin. 2001;17(1):107-124.

Krinsley

JS.

Association between hyperglycemia and increased hospital mortality in a heterogeneous population of critically ill patients. Mayo Clin Proc 2003;78(12):1471-1478.

Finney

Zekveld

Elia

Evans

TW.

Glucose control and mortality in critically ill patients. JAMA. 2003;290(15):2041-2047.

Christiansen

Toft

Jørgensen

Andersen

Tønnesen

Hyperglycaemia and mortality in critically ill patients. A prospective study. Intensive Care Med. 2004;30(8):1685-1688.

Krinsley

JS.

Effect of an intensive glucose management protocol on the mortality of critically ill adult patients. Mayo Clin Proc. 2004;79(8):992-1000.

10.

Chase

Shaw

Le Compte

et al . Implementation and evaluation of the SPRINT protocol for tight glycaemic control in critically ill patients: a clinical practice change. Crit Care. 2008;12(2):R49.

11.

Brunkhorst

Engel

Bloos

et al . Intensive insulin therapy and pentastarch resuscitation in severe sepsis. N Engl J Med. 2008;358(2):125-139.

12.

Hirasawa

Oda

Nakamura

Blood glucose control in patients with severe sepsis and septic shock. World J Gastroenterol. 2009;15(33):4132-4136.

13.

Finfer

Chittock

et al . Intensive versus conventional glucose control in critically ill patients. N Engl J Med. 2009;360(13):1283-1297.

14.

Wiener

Larson

RJ.

Benefits and risks of tight glucose control in critically ill adults: a meta-analysis. JAMA. 2008;300(8):933-944.

15.

Egi

Bellomo

Stachowski

et al . Hypoglycemia and outcome in critically ill patients. Mayo Clin Proc. 2010;85(3):217-224.

16.

Bagshaw

Bellomo

Jacka

et al . The impact of early hypoglycemia and blood glucose variability on outcome in critical illness. Crit Care. 2009;13(3):R91.

17.

NICE-SUGAR Study Investigators. Hypoglycemia and risk of death in critically ill patients. N Engl J Med. 2012;367(12):1108-1118.

18.

Krinsley

Preiser

JC.

Time in blood glucose range 70 to 140 mg/dl >80% is strongly associated with increased survival in non-diabetic critically ill adults. Crit Care. 2015;19:179.

19.

Chase

Andreassen

Jensen

Shaw

GM.

Impact of human factors on clinical protocol performance: a proposed assessment framework and case examples. J Diabetes Sci Technol. 2008;2(3):409-416.

20.

Chase

Le Compte

Suhaimi

Tight glycemic control in critical care—the leading role of insulin sensitivity and patient variability: a review and model-based analysis. Comput Methods Programs Biomed. 2011;102(2):156-171.

21.

Suhaimi

Le Compte

Preiser

et al . What makes tight glycemic control tight? The impact of variability and nutrition in two clinical studies. J Diabetes Sci Technol. 2010;4(2):284-298.

22.

Signal

Pretty

Chase

Le Compte

Shaw

GM.

Continuous glucose monitors and the burden of tight glycemic control in critical care: can they cure the time cost?

J Diabetes Sci Technol. 2010;4(3):625-635.

23.

Pretty

Le Compte

Chase

et al . Variability of insulin sensitivity during the first 4 days of critical illness: implications for tight glycemic control. Ann Intensive Care. 2012;2(1):17.

24.

Holzinger

Warszawska

Kitzberger

et al . Real-time continuous glucose monitoring in critically ill patients: a prospective randomized trial. Diabetes Care. 2010;33(3):467-472.

25.

Boom

Sechterberger

Rijkenberg

et al . Insulin treatment guided by subcutaneous continuous glucose monitoring compared to frequent point-of-care measurement in critically ill patients: a randomized controlled trial. Crit Care. 2014;18(4):453.

26.

Facchinetti

Del Favero

Sparacino

Castle

Ward

Cobelli

Modeling the glucose sensor error. IEEE Trans Biomed Eng. 2014;61(3):620-629.

27.

Reifman

Rajaraman

Gribok

Kenneth Ward

Predictive monitoring for improved management of glucose levels. J Diabetes Sci Technol. 2007;1(4):478-486.

28.

Zimmermann

Lehmann

Hofer

et al . Design of a prospective clinical study on the agreement between the Continuous GlucoseMonitor, a novel device for CONTinuous ASSessment of blood GLUcose levels, and the RAPIDLab (R) 1265 blood gas analyser: The CONTASSGLU study. BMC Anesthesiol. 2012;12:24.

29.

Kuure-Kinsey

Palerm

Bequette

BW.

A dual-rate Kalman filter for continuous glucose monitoring. Conf Proc IEEE Eng Med Biol Soc. 2006;1:63-66.

30.

Breton

Kovatchev

Analysis, modeling, and simulation of the accuracy of continuous glucose sensors. J Diabetes Sci Technol. 2008;2(5):853-862.

31.

Lunn

Wei

Hovorka

Fitting dynamic models with forcing functions: application to continuous glucose monitoring in insulin therapy. Stat Med. 2011;30(18):2234-2250.

32.

Signal

Gottlieb

Le Compte

Chase

JG.

Continuous glucose monitoring and trend accuracy: news about a Trend Compass. J Diabetes Sci Technol. 2014;8(5):986-997.

33.

Preiser

Chase

Hovorka

et al . Glucose control in the ICU: a continuing story. J Diabetes Sci Technol. 2016;10(6):1372-1381.

34.

Kovatchev

Breton

Man

Cobelli

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.

35.

Patek

Bequette

Breton

et al . In silico preclinical trials: methodology and engineering guide to closed-loop control in type 1 diabetes mellitus. J Diabetes Sci Technol. 2009;3(2):269-282.

36.

Facchinetti

Sparacino

Cobelli

Modeling the error of continuous glucose monitoring sensor data: critical aspects discussed through simulation studies. J Diabetes Sci Technol. 2010;4(1):4-14.

37.

Crane

Barwell

Gopal

et al . The development of a continuous intravascular glucose monitoring sensor. J Diabetes Sci Technol. 2015;9(4):751-761.

38.

Burmeister

Arnold

MA.

Accuracy of the YSI Stat Plus Analyzer for glucose and lactate. Analytical Lett. 1995;28(4):581-592.

39.

Erickson

Wilding

Evaluation of a novel point-of-care system, the i-STAT portable clinical analyzer. Clin Chem. 1993;39(2):283-287.

40.

Chase

Suhaimi

Penning

et al . Validation of a model-based virtual trials method for tight glycemic control in intensive care. Biomed Eng Online. 2010;9:84.

41.

Hovorka

Chassin

Ellmerer

Plank

Wilinska

ME.

A simulation model of glucose regulation in the critically ill. Physiol Meas. 2008;29(8):959-978.

42.

Chase

Le Compte

Preiser

Shaw

Penning

Desaive

Physiological modeling, tight glycemic control, and the ICU clinician: what are models and how can they affect practice?

Ann Intensive Care. 2011;1(1):11.

43.

Chase

Desaive

Preiser

J-C

. Virtual patients and virtual cohorts: a new way to think about the design and implementation of personalized ICU treatments. In: Vincent J-L ed. Annual Update in Intensive Care and Emergency Medicine. New York, NY: Springer; 2016:435-448.

44.

Damiano

McKeon

El-Khatib

Zheng

Nathan

Russell

SJ.

A comparative effectiveness analysis of three continuous glucose monitors: the Navigator, G4 Platinum, and Enlite. J Diabetes Sci Technol. 2014;8(4):699-708.

45.

Facchinetti

Sparacino

Guerra

Real-time improvement of continuous glucose monitoring accuracy: the smart sensor concept. Diabetes Care. 2013;36(4):793-800.

46.

Kovatchev

Heinemann

Anderson

Clarke

Comparison of the numerical and clinical accuracy of four continuous glucose monitors. Diabetes Care. 2008;31(6):1160-1164.

47.

Keenan

Cartaya

Mastrototaro

JJ.

Accuracy of a new real-time continuous glucose monitoring algorithm. J Diabetes Sci Technol. 2010;4(1):111-118.

48.

Freckmann

Pleus

Link

Zschornack

Klötzer

Haug

Performance evaluation of three continuous glucose monitoring systems: comparison of six sensors per subject in parallel. J Diabetes Sci Technol. 2013;7(4):842-853.

49.

Geoffrey

Brazg

Richard

FreeStyle navigator continuous glucose monitoring system with TRUstart algorithm, a 1-hour warm-up time. J Diabetes Sci Technol. 2011;5(1):99-106.

50.

Stewart

Thomas

Pretty

Chase

Shaw

How should we interpret retrospective blood glucose measurements? Sampling and interpolation. Paper presented at: 20th World Congress of the International Federation of Automatic Control; Toulouse, France, 2017.

51.

Blixt

Rooyackers

Isaksson

Wernerman

Continuous on-line glucose measurement by microdialysis in a central vein. A pilot study. Crit Care. 2013;17(3):R87.

52.

Luijf

Mader

Doll

et al . Accuracy and reliability of continuous glucose monitoring systems: a head-to-head comparison. Diabetes Technol Ther. 2013;15(8):722-727.

53.

Zisser

Bailey

Schwartz

Ratner

Wise

Accuracy of the SEVEN continuous glucose monitoring system: comparison with frequently sampled venous glucose measurements. J Diabetes Sci Technol. 2009;3(5):1146-1154.