Fault and meal detection by redundant continuous glucose monitors and the unscented Kalman filter

Biomedical Signal Processing and Control 38 (2017) 86–99

Contents lists available at ScienceDirect

Biomedical Signal Processing and Control journal homepage: www.elsevier.com/locate/bspc

Fault and meal detection by redundant continuous glucose monitors and the unscented Kalman filter Zeinab Mahmoudi a,b , Kirsten Nørgaard c , Niels Kjølstad Poulsen a , Henrik Madsen a , John Bagterp Jørgensen a,∗ a

Department of Applied Mathematics and Computer Science, Technical University of Denmark, DK-2800 Kgs. Lyngby, Denmark Danish Diabetes Academy, Odense University Hospital, DK-5000 Odense, Denmark c Department of Endocrinology, Copenhagen University Hospital Hvidovre, DK-2650 Hvidovre, Denmark b

a r t i c l e

i n f o

Article history: Received 6 December 2016 Received in revised form 7 April 2017 Accepted 8 April 2017 Keywords: Continuous glucose monitoring sensor Fault detection Unscented Kalman filter Adaptive filtering Sensor redundancy

a b s t r a c t The purpose of this study is to develop a method for detecting and compensating the anomalies of continuous glucose monitoring (CGM) sensors as well as detecting unannounced meals. Both features, sensor fault detection/correction and meal detection, are necessary to have a reliable artificial pancreas. The aim is to investigate the best detection results achievable with the proposed detection configuration in a perfect situation, and to have the results as a benchmark against which the imperfect scenarios of the proposed fault detection can be compared. The perfect situation that we set up here is in terms of a patient simulation model, where the model in the detector is the same as the patient simulation model used for evaluation of the detector. The detection module consists of two CGM sensors, two fault detectors, a fault isolator, and an adaptive unscented Kalman filter (UKF). Two types of sensor faults, i.e., drift and pressure induced sensor attenuation (PISA), are simulated by a Gaussian random walk model. Each of the fault detectors has a local UKF that receives the signal from the associated sensor, detects faults, and finally tunes the adaptive UKF. A fault isolator that accepts data from the two fault detectors differentiates between a sensor fault and an unannounced meal appearing as an anomaly in the CGM data. If the fault isolator indicates a sensor fault, a method based on the covariance matching technique tunes the covariance of the measurement noise associated with the faulty sensor. The main UKF uses the tuned noise covariances and fuses the CGM data from the two sensors. The drift detection sensitivity and specificity are 80.9% and 92.6%, respectively. The sensitivity and specificity of PISA detection are 78.1% and 82.7%, respectively. The fault detectors can detect 100 out of 100 simulated drifts and 485 out of 500 simulated PISA events. Compared to a nonadaptive UKF, the adaptive UKF reduces the deviation of the CGM measurements from their paired blood glucose concentrations from 72.0% to 12.5% when CGM is corrupted by drift, and from 10.7% to 6.8% when CGM is corrupted by PISA. The fault isolator can detect 199 out of 200 unannounced meals. The average change in the glucose concentrations between the meals and the detection time points is 46.3 mg/dL. © 2017 Elsevier Ltd. All rights reserved.

1. Introduction As the continuous glucose monitoring (CGM) sensor is one of the main parts of an artificial pancreas (AP), the anomalies and faults associated with the CGM sensor affect the performance and reliability of the AP considerably. Consequently, the safety of patients in an AP is critically connected to the reliability of the CGM sensors

∗ Corresponding author. E-mail addresses: [email protected] (Z. Mahmoudi), [email protected] (K. Nørgaard), [email protected] (N.K. Poulsen), [email protected] (H. Madsen), [email protected] (J.B. Jørgensen). http://dx.doi.org/10.1016/j.bspc.2017.05.004 1746-8094/© 2017 Elsevier Ltd. All rights reserved.

[1–4]. Recent advances in sensor technology have made it possible to achieve high accuracy CGM sensors. An example is the CGM sensor under development by Roche (Roche Diagnostics GmbH, Mannheim, Germany), which provides an aggregated mean absolute relative difference (MARD) between the CGM readings and their paired capillary blood glucose (BG) measurements of 9.2% [5]. A study by Bailey et al. [6] also indicated that the FreeStyle Libre Flash glucose monitoring system (Abbott Diabetes Care, Alameda, CA) has an overall MARD of 11.4% without the need for calibration within two weeks. However, the Flash sensor is not suited for the AP application. Dexcom (San Diego, CA) also tested the ® accuracy of the G4 Platinum (G4P) CGM system with the 505

Z. Mahmoudi et al. / Biomedical Signal Processing and Control 38 (2017) 86–99

algorithm, achieving a MARD of 13% [7]. Furthermore, the Dexcom ® G5 mobile CGM has been recently approved for insulin-dosing in ® Europe, which indicates the acceptable reliability of the G5 system for nonadjunctive use [8]. Despite these advances in the sensor technology and accuracy, the CGM fault detection is still a challenge that requires attention and it serves as an active research area [9]. In this study, we enhance the CGM data by developing an algorithm for fault and meal detection. Fault detection improves the CGM safety and detection of unannounced meals helps developing a fully automated AP. The faults that we aim to detect are drift and the pressure induced sensor attenuation (PISA) artifact. Drift remains the main reason for sensor calibration [10], and calibration for at least twice a day is recommended to compensate for the sensor drift. Drift is the falsely slow variations of the CGM readings due to foreign body response causing the inflammatory cells to migrate to the sensor insertion site [11]. The inflammatory cells produce compounds that interact with the glucose sensor and erroneously increase or decrease the glucose readings over the course of time. PISA is another sensor anomaly and it is the low signal readings caused by the compression of the sensor insertion site due to pressure on the sensor [12]. Several methods have been proposed for outlier detection and smoothing the CGM data. These methods include finite and infinite impulse response filters [13,14], as well as model-based outlier detection with adaptive and nonadaptive Kalman filters [15–18]. While these denoising methods effectively smooth the CGM signal by filtering out the spiky outliers and noise, they are insufficient to remove the manifestations of the sensor artifacts, such as drift, from the CGM data. The reason is that the time scale of a sensor drift could be comparable with the time scale of the signature of the physiological events such as meal, insulin, exercise, and stress on the CGM signal. This makes it difficult to detect the sensor drift and differentiate it from the metabolic changes in the BG. Consequently, a robust and safe AP requires additional sources of information through physical redundancy (using more than one CGM sensor) and analytical redundancy (using a model in addition to the measurements) to detect and discriminate between drift and the physiological events such as unannounced meals. While previous studies suggest methods for drift and PISA detection [19,20], the literature on filtering out the drift and PISA events upon detection is sparse. We previously developed a fault detector for the CGM data using the extended Kalman filter (EKF) [21], which we further improved by using an unscented Kalman filter (UKF) [22]. In the current study, we use the fault detector with the UKF for the detection of drift, PISA, and unannounced meals, and also for tuning an adaptive UKF upon fault detection, in a two-CGM sensor configuration. Fig. 1 depicts the schematic configuration of the two CGM sensors, the insulin pump, and the signal processing unit, for the proposed method in a single-hormone AP. The adaptive UKF allows modifying the covariance of the measurement noise in case the CGM sensors are fault-corrupted. This mitigates the signature of a fault on the filtered CGM data. We also designed a fault isolator unit to discriminate between a sensor fault and a change in the CGM signal due to the patient’s metabolism variations. Examples of metabolic variations are meal, exercise, stress, and variation of the insulin sensitivity. Any of these events, if not correctly announced to the model, manifests itself in the CGM signal in a way that is similar to the manifestation of a sensor anomaly. As different treatment should be applied for a sensor anomaly and a metabolic change, it is also important that these two categories of events are detected and differentiated accordingly. In our method, the UKF is tuned only if the fault isolator indicates a sensor fault. The rest of the paper is structured as follows. Section 2 describes the Medtronic virtual patient (MVP) model for type 1 diabetes [23], and our methods for simulating sensor drift and PISA. The MVP

87

Insulin pump

Plaorm: - Control algorithm - CGM monitoring and fault detecon - Data integraon

CGM sensor 1 CGM sensor 2 Fig. 1. Configuration of the two continuous glucose monitoring sensors, the insulin pump, and the signal processing platform-implemented in a smart phone – for an artificial pancreas.

CHO absorpon

PD d 1

2

Insulin absorpon

PK

Fig. 2. The Medtronic virtual patient model for type 1 diabetes [23]. d: CHO ingestion rate; d1 and d2 : glucose masses (mg) in the stomach and small intestine compartments, respectively; RA : glucose appearance rate; GB : blood glucose concentration; GI : interstitial glucose concentration; u: insulin input rate; ISC : SC insulin concentration; Ip : plasma insulin concentration; Ieff : insulin effect; : measurement noise; y: CGM measurement.

model is used for patient simulation and also in the UKFs. Section 3 explains the UKF for filtering, fault detection, and fault isolation and meal detection. At the end, Section 3 presents the method for tuning the measurement noise covariance with the local UKFs of the fault detectors. Section 4 presents the results followed by the discussions in Section 5 and concluding remarks in Section 6. 2. Simulation 2.1. Virtual patient model We simulated the MVP model [23], and the two-compartmental model of carbohydrate (CHO) absorption [24], by means of stochastic differential equations (SDEs). The model is based on the work by Hovorka and associates [24–26]. Fig. 2 depicts the model and Appendix A describes the model equations. In Fig. 2, u is the subcutaneous (SC) insulin input rate (␮U/min) and contains both basal and bolus insulin administrations. ISC , Ip , and Ieff are the SC insulin concentration (mU/L), the plasma insulin concentration (mU/L), and the effect of insulin (min−1 ), respectively. GB is the blood glucose concentration. Glucose diffuses from capillary blood into interstitial fluid and the CGM sensor measures the interstitial glucose (GI ) concentration (mg/dL). The input d denotes the CHO ingestion rate (g/min), d1 and d2 are the glucose masses (mg) in the stomach and small intestine compartments, and RA is the glu-

88

Z. Mahmoudi et al. / Biomedical Signal Processing and Control 38 (2017) 86–99

cose appearance rate (mg/dL/min).  is the measurement noise that is defined in the next section.

400 350 C1

CGM measurement (mg/dL)

2.2. Fault simulation Eq. (1) presents a model for simulating CGM data. According to this model, the CGM measurement, y, consists of the interstitial glucose concentration, GI , and an additive discrete-time measurement noise, : yk = GI,k + k ,

(1)

where k denotes the time index. When the sensor is fault-free,  = n , in which n is the nominal measurement noise with covariance Rn . The noise n , identified by Facchinetti et al. [27], consists of two autoregressive processes according to

300

A

E1 B

250 B

200 180 D1

D2

D3

100 70 E2

ˆ k, n,k = ck + ϑ

(2a)

ck = 1.23ck−1 − 0.3995ck−2 + ıc,k ,

(2b)

ˆ k = 1.013ϑ ˆ k−1 − 0.2135ϑ ˆ k−2 + ı , ϑ ˆ ϑ,k

(2c)

2

ıc,k ∼N(0, 11.3 mg /dL ), 2 2 ıϑ,k ˆ ∼N(0, 14.45 mg /dL ).

Facchinetti et al. [27] identified  by a CGM model that contains the GB − GI physiology, the sensor calibration, and a zero-mean autoregressive random noise generated by the sensor. The process, c, represents the CGM inaccuracy (with respect to GB ) due to the imperfect model of GB − GI dynamics, and the calibration error. The calibration error is caused by the erroneous BG measurements, used as reference, and also by the suboptimal calibration function. The ˆ is the estimated zero-mean random noise of the sensor. process ϑ When y is corrupted by a fault,  = f , in which f is the faulty measurement noise. It contains both the nominal noise and the fault. The noise, f , is defined by (3)

where f is the additive sensor fault, i.e., drift or PISA. A faulty CGM reading can be also presented in the form yf,k = yn,k + fk ,

70

130

180

240

290

400

Blood glucose concentration (mg/dL) Fig. 3. Modified Clarke error grid classification zones, used for simulating the sensor drift. Zone A contains the most accurate measurements. The CGM accuracy degrades as the measurements fall in zones B, C, D, and E. The CGM readings in zone E are the least accurate.

2

f,k = n,k + fk ,

C2

(4)

in which yf is the sensor measurement with the faulty noise, and yn is the CGM measurement with the nominal noise.

Table 1 Sub-zone definition for simulating anomaly F. Sub-zone

F∼

C1 : 70 < GB ≤ 290 C2 : 130 < GB ≤ 180 D1 : 0 < GB ≤ 175/3 D2 : 175/3 < GB ≤ 70 D3 : 240 < GB ≤ Ul E1 : 0 < GB ≤ 70 E2 : 180 < GB ≤ Ul

U((GB + 110 − yn ), (Ul − yn )] U(−yn , (7/5 · GB − 182 − yn )] U((70 − yn ), (180 − yn )] U((6/5 · GB − yn ), (180 − yn )] U((70 − yn ), (180 − yn )] U((180 − yn ), (Ul − yn )] U(−yn , (70 − yn )]

Ul : the upper bound of the glucose concentration that sensor can measure.

where s is the onset of the drift, and m is the drift duration, which is 12 h. The drift ends at time sample s + m, indicating sensor calibration. The onset of the drift is chosen randomly. A Gaussian random walk models the drift, specified by dk = dk−1 + d. d∼N(, d2 ), ds = 0.

(6)

The need for calibration every 12 h implies that the mathematical expectation of d approaches the anomalous value F at sample s + m:

2.2.1. Drift For simulating drift in the CGM data, we used the definition of the CGM clinical accuracy provided by the Clark Error Grid Analysis (CEGA). CEGA is an accuracy assessment tool that quantifies the accuracy of the CGM readings using the BG or the plasma glucose (PG) concentrations as reference. CEGA quantifies the risk of a CGM-based decision making for the patient according to five classification zones. The CGM readings falling in zone A and zone B are clinically acceptable, with the CGM measurements in zone A being the most accurate readings. The CGM readings in zones C, D and E are not reliable. Drift is an additive value that gradually brings y to any of the zones C, D, and E. Fig. 3 illustrates our modification to the CEGA by dividing zone C into C1 and C2 , zone D into D1 , D2 and D3 , and zone E into E1 and E2 . Let d denote the sensor drift and let yd denote the driftcorrupted CGM measurement. The summation of the nominal CGM measurement and the simulated drift generates the drifty CGM measurement, defined as

When F is added to yn , the measurement falls in one of the zones C, D or E. Table 1 gives the sub-zone definitions and the uniform distributions for simulating F. The value of F depends on the blood glucose concentration, GB , and the measurement, yn , at sample s + m, according to the uniform distributions in the table. When more than one sub-zone is possible for F, a random permutation chooses the sub-zone from the possible options. d2 in (6) is set to 1 mg2 /dL2 . Because d has a normal distribution, the following holds:

yd,k = yn,k + dk , k = s + 1, s + 2, . . ., s + m,

=

(5)

E{ds+m } = F.

E{ds+m } = m.

(7)

(8)

Eqs. (7) and (8) yield  according to F . m

(9)

Z. Mahmoudi et al. / Biomedical Signal Processing and Control 38 (2017) 86–99

89

Nonlinear Transformaon

Fig. 5. The principle of the unscented transformation. The sigma points approximate the probability distribution and the nonlinear function transforms the sigma points. The statistics of the transformed sigma points approximate the nonlinearly transformed distribution.

Fig. 4. PISA in the real data confirmed by the SMBG measurement. The data are col® lected by the Medtronic Paradigm VeoTM system. The data are recorded in mmol/L. They are converted to mg/dL by multiplication by 18, i.e. 3.3 mmol/L corresponds to 60 mg/dL, 3.9 mmol/L corresponds to 70 mg/dL, 5 mmol/L corresponds to 90 mg/dL, and 10 mmol/L corresponds to 180 mg/dL.

2.2.2. Pressure induced sensor attenuation To simulate PISA, we use PISA information from the outpatient CGM data of 10 patients with type 1 diabetes wearing the Medtronic ® Paradigm VeoTM CGM augmented insulin pump (Medtronic Minimed, Inc., Northridge, CA). The data were downloaded for regular screening in Copenhagen University Hospital Hvidovre. The average number of reviewed data is 72 days per patient. The CGM data are accompanied by insulin and self-monitored BG (SMBG) information. With the supervision of an experienced diabetologist, we inspected the data to detect PISA events. A PISA was detected if it was confirmed by a SMBG measurement, i.e., if the nadir of the PISA event was considerably lower than its concurrent SMBG value. This approach confirmed 18 PISA events. The small number of detected PISA events is due to the rare confirmatory SMBG measurements. Fig. 4 illustrates an example of a confirmed PISA event in the actual data of a patient with type 1 diabetes. The duration of a PISA event is the time interval between the onset of the signal fall (when the slope starts to be negative) and the time point that the signal returns from the nadir to a non-increasing state (when the slope becomes zero). The PISA nadir is the difference between the minimum value of the PISA event and the CGM value immediately before the onset of the fall. The nadir is located centrally. The mean (SD) of the detected PISA duration is 51.1 (10.8) min, and the mean (SD) of the nadir amplitude is −116.2 (41.4) mg/dL. Sampling from the normal distribution of N(51.1 mg/dL, 10.82 mg2 /dL2 ) gives the duration and sampling from the normal distribution of N(−116.2 mg/dL, 41.42 mg2 /dL2 ) yields the nadir of a simulated PISA event. Having the duration and nadir of the event, the following procedure simulates a PISA. A Gaussian random walk with the step size distribution of N(1 mg/dL, 1 mg2 /dL2 ), where 1 = nadir/(0.5 duration), models the PISA event from the onset to nadir. To complete the event from the nadir to the end, a Gaussian random walk with the step size distribution of N(2 mg/dL, 1 mg2 /dL2 ), where 2 =− 1 , is used. 3. The unscented Kalman filter

(sigma points) are chosen to represent the mean and covariance and potentially higher order moments of the states [28–30]. The nonlinear function f is applied to the sigma points and the statistics of the transformed points form an estimate of the nonlinearly transformed distribution. 3.1. Model formulation The MVP model for the patients with type 1 diabetes in [23] can be extended with the diffusion term, and it can be presented in the SDE form dx(t) = f (x(t), u(t), d(t))dt +  · dω(t),

(10)

where f is the nonlinear state-space model, x is the vector of the state variables, i.e., x = [ISC Ip Ieff GB GI d1 d2 ]T , u is the SC insulin input rate, d is the CHO ingestion rate,  is the diffusion coefficient, and ω is a standard Wiener process. Appendix A describes the equations of the model. The sensor measurements can be modeled as yk = g(xk ) + k ,

(11)

where y is the CGM measurement, and  is a zero-mean Gaussiandistributed discrete-time measurement noise with covariance R. 3.2. Multi-step prediction (time update) This section presents the multi-step prediction with the UKF. The prediction steps are j = 1, 2, . . . l, where l is the prediction horizon. The associated weights of the 2n + 1 sigma points, in which n is the dimension of the state space, are calculated as , (n + )

(12a)

, {(n + )(1 − ˛2 + ˇ)}

(12b)

(i)

1 , i = 1, . . ., 2n, {2(n + )}

(12c)

(i)

1 , i = 1, . . ., 2n, {2(n + )}

(12d)

(0)

Wm = (0)

Wc

=

Wm = Wc =

(0)

(2n) T

wm = [Wm . . .Wm ] , wc =

(12e)

(0) (2n) [Wc . . .Wc ],

(12f) T

We used the UKF for the fault detection and adaptive filtering of the CGM measurements. The intuition behind the unscented transformation in the UKF is that it is easier to approximate a probability distribution than it is to approximate a nonlinear function. Fig. 5 illustrates the unscented transformation. A set of points

W = (I − [wm . . .wm ]) × diag(wc ) × (I − [wm . . .wm ]) ,

(12g)

where = ˛2 (n + ) − n is a scaling parameter, the parameter ˛ ∈ (0, 1] influences how far the sigma points are spread around the mean (here we set ˛ = 0.01), and ˇ is a positive constant (ˇ = 2 is optimal for the Gaussian distribution). The parameter is often set to 0.

90

Z. Mahmoudi et al. / Biomedical Signal Processing and Control 38 (2017) 86–99

ˆ are generated via a deterministic sampling The sigma points, X, of the distribution using the mean, xˆ , and the covariance, P, according to the following: Xˆ k+j−1 =

[ˆxk+j−1|k . . . xˆ k+j−1|k ]   √   + c 0 Pk+j−1|k − Pk+j−1|k (0)

(2n)

= [mx

. . . mx

dt

T Pk+1|k+1 = Pk+1|k − Kk+1 Sk+1 Kk+1 .



t ∈ [tk+j−1 tk+j ],

(13b)

Xˆ k+j = Xˆ k+j−1 (tk+j ).

(13c)

The weighted average of the transformed sigma points gives the mean of the predicted state according to xˆ k+j|k = Xˆ k+j wm .

(13d)

Integrating the following differential equation yields the propagated covariance, Pk+j|k : dt

=





Xˆ k+j−1 (t)Wf T Xˆ k+j−1 (t), u(t), d(t)





T +f Xˆ k+j−1 (t), u(t), d(t) W Xˆ k+j−1 (t) +  T ,

Pk+j|k = Pk+j−1 (tk+j ).

(13f)

The measurement model transforms the new sigma points: (0) (2n) Yˆ k+j = g(X˜ k+j ) = [my . . . my ].

(13h)

The weighted average of the measurement sigma points gives the predicted measurement as yˆ k+j|k = Yˆ k+j wm .

(13i)

The prediction error is ek+j|k = yk+j − yˆ k+j|k ,

(13j)

and the covariance of the j-step prediction error is (13k)

The cross-covariance matrix of the state and measurement is determined by xy

(13l)

3.3. Filtering (measurement update) The filtering step in the UKF is similar to the linear Kalman filter equations. The filter gain is computed as Kk+1 =

xy Sk+1

−1 Sk+1 ,

Fig. 6 illustrates an overview of the proposed method in the study. The two UKFs of the fault detectors are local state estimators with the fixed process and measurement noise covariances. The purpose of each filter is to detect a fault in one of the sensors. The faults being detected are drift and PISA. The fault isolator differentiates between a sensor fault and a change in the patient’s metabolism. We considered a change in the metabolism to be an unannounced meal. If the detected anomaly is not a sensor fault, the adaptive UKF uses the nominal measurement noise covariances. The adaptive UKF takes measurements from both sensors. The measurement noise covariances in the adaptive UKF have either the nominal values or the tuned values, depending on the fault isolator outcome. The process noise covariance remains unchanged, and it is unaffected by the state of fault and meal detection.

3.5. Fault detection

To reduce the prediction bias, a new set of sigma points is generated from the predicted mean and covariance of the state according to   √ X˜ k+j = [ˆxk+j|k . . . xˆ k+j|k ] + c [0 Pk+j|k − Pk+j|k ] (13g) (0) (2n) ˜ x ... m ˜ x ]. = [m

T Sk+j = X˜ k+j W Yˆ k+j .

Appendix B outlines the filter parameters Q and R for the fault detectors and the main UKF.

(13e)

t ∈ [tk+j−1 tk+j ],

T Sk+j = Yˆ k+j W Yˆ k+j + R.

(14c)

3.4. General description of the method

(t) = f Xˆ k+j−1 (t), u(t), d(t) ,

dPk+j−1 (t)

(14b)

Finally, the updated covariance of the state is given by

],



xˆ k+1|k+1 = xˆ k+1|k + Kk+1 ek+1|k .

(13a)

where c = ˛2 (n + ). Because Pk+j−1|k is positive definite, we use a Cholesky decomposition to obtain the matrix square root, and  Pk+j−1|k is the lower triangular matrix of chol(Pk+j−1|k ). The nonlinear function propagates each of the sigma points: dXˆ k+j−1

and the state mean, updated by the measurement, is

(14a)

The essence of the presented fault detection is the multi-step prediction in a window of size l, by means of the UKF, and applying a statistical test on the prediction residuals. For the UKFs in the fault detectors, the noise  in the measurement model (11) is set to n in (2). This implies that in case of a faulty y signal, the model of the measurement noise does not fit the observed measurement noise. The prediction error manifests this mismatch which is used for fault detection. denote the CGM measurements and l denote the preLet {yk }N k=1 diction horizon. Let the time indices be k = 0, 1, . . ., N − l, and the prediction index be 1 ≤ j ≤ l. This implies that 1 ≤ k + j ≤ N. Let k+l denote the vector of the prediction errors in the l-sample fault detection window, defined by

⎡

ek+1|k

⎤

⎡

yk+1 − yˆ k+1|k

⎤

⎢e ⎥ ⎢ ⎥ ⎢ k+2|k ⎥ ⎢ yk+2 − yˆ k+2|k ⎥ ⎥=⎢ ⎥. ⎥ ⎢ ⎥ .. ⎣ ... ⎦ ⎣ ⎦ .

k+l = ⎢ ⎢

ek+l|k

(15)

yk+l − yˆ k+l|k

In every prediction window, the local UKF in each of the fault detectors generates the vector k+l from the multi-step predictions and the sensor measurements. Fig. 7 illustrates this process. Analogously to a linear system, the cross-covariances of the multi-step prediction errors may be computed as [31]

Sk+i,k+j = ek+i|k , ek+j|k  =

⎧ Yˆ W Yˆ T , ⎪ ⎪ ⎨ k+i k+j

T +R , Yˆ k+i W Yˆ k+j n

⎪ ⎪ ⎩ Yˆ

ˆT k+j W Yk+i ,

i > j, i = j, i < j.

(16)

Z. Mahmoudi et al. / Biomedical Signal Processing and Control 38 (2017) 86–99

91

Fig. 6. Detection of sensor faults and unannounced meals by means of the UKF and sensor redundancy.

Fig. 7. Construction of the prediction error vector, k+l , from the multi-step predictions and the CGM measurements in the l-sample fault detection window.

Under the fault-free condition, the relation between T2 and the F-distribution is determined by [33]

The covariance matrix of k+l is calculated as

⎡

Sk+1,k+1

Sk+1,k+2

···

Sk+1,k+l

⎤

⎢S ⎥ ⎢ k+2,k+1 Sk+2,k+2 · · · Sk+2,k+l ⎥ ⎢ ⎥. k+l = k+l , k+l  = ⎢ ⎥ .. .. .. .. ⎣ ⎦ . . . . Sk+l,k+1

Sk+l,k+2

···

T 2∼ (17)

Sk+l,k+l

In the absence of fault and under the assumption that the measurement noise  is Gaussian, k+l has the multivariate normal distribution of Nl ([0, 0, . . . 0]T , k+l ). This is the null hypothesis (H0 ) of our fault detection. To test the null hypothesis, in every l-sample fault detection window, we form the Hotelling’s T2 test statistic as [32] T 2 = Tk+l −1  . k+l k+l

(18)

l(q2 − 1) F , q(q − l) l,q−l

(19)

where Fl,q−l is the F-distribution with degrees of freedom l and q − l, and q is the number of measurement data points used for estimation of k+l . If q is sufficiently large, the above distribution can be approximated by a 2 -distribution with l degrees of freedom, and we have [34] T 2 ∼ 2 (l).

(20)

To verify whether q is large enough, we refer to (17) which indicates that the entries of the matrix are the estimates of S. According to (16), Yˆ contributes to the estimation of S, and back˜ The sigma ing to (13h), Yˆ is computed based on the estimate of X. points X˜ are themselves formed by the xˆ k+j|k predictions according to (13g), and xˆ k+j|k is estimated based on the measurement data points up to k, namely {yi }ki=1 . Therefore q = k, and if k »1, the con-

92

Z. Mahmoudi et al. / Biomedical Signal Processing and Control 38 (2017) 86–99

dition for using (20) is fulfilled. We assumed the condition k »1 is satisfied when k > 3l. Based on this assumption, we made the following decision logic for the fault detection:



H0 is retained (fault-free CGM),

if T 2 ≤ 2 (˛, l),

H0 is rejected (faulty CGM),

if T 2 > 2 (˛, l),

(21)

where 2 (˛, l) is the critical value of the 2 -distribution with l degrees of freedom at the significance level of ˛ (˛ = 0.05 is used here). The fault detector announces a fault if the null hypothesis is rejected. When a fault is detected, a flag is lifted for the interval k + l to k + l + w, where w is the length of the moving step of the sliding fault detection window. The next detection window starts at sample k + w, indicating l − w samples overlap with the last window. Each of the fault detectors in Fig. 6 follows the described method. For the drift detection, both fault detectors have l = 90 min and w = 15 min, and for the PISA detection they have l = 10 min and w = 1 min. Fig. 8 indicates the fault detection method.

3.6. Detection of fault sign In the l-sample fault detection window, if yk+j is larger than the upper limit of the 95% CI of yˆ k+j|k prediction, a positive sign is assigned to the fault at sample k + j. Similarly, if yk+j is smaller than the lower limit of the 95% CI of yˆ k+j|k prediction, the fault at sample k + j is assigned a negative sign. The fault isolator incorporates the detected sign, as described in the next section.

Table 2 Drift detection for 50 days of simulated data for each sensor. Total number of drifts Number of detected drifts (true positive) Sensitivity% = TP duration/drift duration Specificity% = 1 − FP duration/data duration without drift Detection delay (min) = delay between the start of the drift and the onset of the drift detection flag Sensitivity of sign detection% BG-CGM deviation at the detection time = (|GB,k − yk |/GB,k ) %, k = detection time a

100% 38.2 (7.7)

3.8. Tuning of the adaptive UKF If the fault isolator indicates that one or the two of the sensors are faulty, an adaptive law, based on the covariance matching technique, tunes the covariance (s) of the measurement noise for the faulty sensor (s). Maybeck [35,Chapter 10] presents the covariance matching technique to derive the adaptive law for the noise covariances. We use this technique in the l-step prediction mode, according to (22). yy,(i)

Sk+l







(i) (i) = Yˆ k+l W Yˆ k+l yy,(i)

where Sk+l

T (22a)

,

is covariance of the l-step predicted y for sensor i.

(i) Sk+l

(i)

is the theoretical covariance of ek+l|k for sensor i, and it is defined as (i)

The basic assumption for fault isolation is that it is very unlikely that the two CGM sensors inserted on two different sites of the body have the same fault magnitude at the same time. Therefore, if both of the CGM signals are announced anomalous by the fault detectors and have the same anomaly magnitude, the anomaly should be attributed to a metabolic change and not a sensor fault. The below mentioned rules implement this idea. Rule 1: if only one of the fault detectors indicates a fault in the signal, the detected anomaly is a sensor fault and the covariance of the measurement noise associated to the faulty sensor is tuned according to the method described later. Rule 2: if both fault detectors indicate a fault, a further test should be applied on the CGM data to isolate the fault. If the number of times that the faults of both sensors have the same sign in the l-sample fault detection window exceeds the threshold of 15, the detected anomaly is a metabolic change, i.e., an unannounced meal, otherwise both sensors are announced faulty. Counting the number of times that the CGM measurements, y, are outside the 95% CI of yˆ k+j|k predictions within the fault detection window gives a rough approximation of the fault magnitude, because it is an indicator of the deviation between the interstitial glucose measurements made by the sensor and the interstitial glucose concentration predicted by the model. Therefore, counting the number of times that the faults of the two sensors have the same sign is a tool to compare the fault magnitude of the two sensors. The specifications of the fault isolator are set to the drift detection mode, as PISA is a distinct sensor fault and does not need to be isolated from a metabolic change. Although a hypoglycemic event induced by the insulin doses much larger than the daily doses may have a similar manifestation to a PISA event, it is rare that patients take such large doses of insulin in the daily life.

94.8 (73.6)

Mean (SD) across the simulation days.

yy,(i)

Sk+l = Sk+l 3.7. Fault isolation and meal detection

73 positive drifts, 27 negative drifts 100 80.9 (13.7)a 92.6 (7.1)

+ R(i) .

(22b) (i)

The sample covariance of ek+l|k is computed as (i) Sˆ k+l =

1 p−1



k+l 

(i)

 

er|r−l

(i)

er|r−l

T ,

(22c)

r=k+l−p+1

where p is the length of the data sequence used for estimating the sample covariance. The length p is set to 30 samples. The tuned measurement covariance Rˆ (i) for sensor i is yy,(i) (i) Rˆ (i) = Sˆ k+l − Sk+l ,

(22d)

which is used only in the adaptive UKF and not in the local UKFs. (i)

The sequence {er|r−l }

k+l r=k+l−p+1

is smoothed with a Savitzky-Golay

filter (SGF) to reduce random noise before covariance estimation [36]. The estimated covariance Rˆ (i) is used in the adaptive UKF for the next w samples corresponding to the moving step of the fault detection sliding window. 4. Results The results are obtained over 50 simulations of one-day data with different realizations of the process and measurement noises. The CGM data of the two sensors have the same realization of the process noise but two different realizations of the nominal measurement noise. The realizations of the sensor faults are also different for the two CGM sensors. Tables 2 and 3 present the results of drift detection and PISA detection, respectively. The true positive (TP) duration is the interval that the fault flag is on, while the fault presents in the data. The interval that the fault flag is on without the presence of fault is considered false positive (FP) duration. The definitions of the performance metrics are the same in Tables 2 and 3.

Z. Mahmoudi et al. / Biomedical Signal Processing and Control 38 (2017) 86–99

93

Fig. 8. The fault detection algorithm, based on the UKF.

Table 3 PISA detection results for 50 days of simulated data for each sensor.a Total number of PISA Duration of the simulated PISA events Depth of the simulated PISA events Number of detected PISA (true positive) Sensitivity% Specificity% Detection delay (min) BG-CGM deviation at the detection time% a

500 51 (10.8) −116.2 (41.4) 485 78.1 (9.6)a 82.7 (3.3) 8.6 (3.2) 25.6 (4.8)

Mean (SD) across the simulation days.

Table 4 Fault isolation (meal detection) false positives. Drift duration (min/day) Mean of type 1 false positive duration (min/day) Mean of type 2 false positive duration (min/day)

720 42.5 19.5

Table 5 Meal detection for 50 days of simulated data and four unannounced meals per day. Number of unannounced meals Number of detected meals TP durationa Change in glucose (mg/dL)b Detection time (min)c

200 199 107.1 (7.2)d 46.3 (21.2) 58.4 (18.7)

a TP is the time duration in which the meal detection flag is on, in the interval between the meal ingestion to 180 min after the meal. b Change in glucose = CGM at detection time − CGM at meal time. c Detection time = The time interval between the meal ingestion and detection point. d Mean (SD) across the simulation days.

Figs. 9 and 10 demonstrate the drift and PISA detections for an example of a one-day simulation. When a calibration is announced, the covariance R is set to the nominal value for the two sensors. An abrupt end in the drift simulates the effect of calibration on the CGM data. Table 4 indicates the false positive information of the fault isolator, which also implies the false positive information for the meal detection. We defined two types of false positive for the fault isolator. False positive type 1 is the time duration when one of the two sensors is drifty and no unannounced meal is present, but the fault isolator indicates an unannounced meal. False positive type 2 is when there is no sensor fault nor any unannounced meal, but the fault isolator indicates the presence of an unannounced meal. Fig. 11 illustrates the performance of the fault isolator for detection of the unannounced meals for four meals per day. To simulate an unannounced meal, we used the meal with the correct content of CHO for the CGM simulation, but we provided the models in the local UKFs and the main UKF with the zero CHO content. Table 5 summarizes the results of applying the meal detection algorithm on 200 unannounced meals.

Table 6 Deviation of the filtered CGM from the simulated BG for the adaptive and nonadaptive UKFs.

Nonadaptive UKF Adaptive UKF a

Fault in the data is drift

Fault in the data is PISA

72.0% (38.7%)a 10.7% (1.3%)

12.5% (3.8%) 6.8% (0.7%)

Mean (SD) across the simulation days.

For comparison purposes and to investigate the degree of CGM accuracy improvement obtained by an adaptive UKF, a nonadaptive UKF, located in the position of the main UKF, filters the CGM in parallel to the adaptive UKF. Table 6 compares the accuracy of the CGM signal, filtered by the adaptive UKF and the nonadaptive UKF, in terms of the MARD of the filtered CGM from the simulated BG. It indicates that using the adaptive UKF results in approximately 50% reduction in the MARD when the CGM signal is PISA-corrupted, and it brings about 85% MARD reduction when the signal is driftcorrupted. 5. Discussions The results indicate that the proposed method can detect 100% of all the simulated drift events, 97% of the simulated PISA events, and almost 100% of all the unannounced meals. The combination of the fault detection and the adaptive filtering decreases the deviation of the filtered CGM from the BG concentrations, as Figs. 9, 10, and Table 6 demonstrate. This improvement in accuracy reduces the need for sensor calibration, because the drift is partially compensated for before the CGM data reach the point that requires sensor calibration. The fault detectors in Fig. 6 are tunable to detect other types of sensor faults. For example, by setting the length of the fault detection window l to 1 in (15) and performing a one-step prediction, the fault detector can detect spiky outliers as we have demonstrated previously [21]. The sensor redundancy provides simultaneous fault and meal detections with high separability between drifts and unannounced meals. This is indicated by the low false positive as Table 4 demonstrates. Turksoy et al. [19] detected CGM drift using principal component analysis based on the large sets of fault-free data used as reference for statistical comparison. By adjusting the fault to signal ratio, they achieved different detection delays and sensitivities. When applied on the simulated data, their method had the detection delay in the range of 10–20 min, which is smaller than the detection delay in our algorithm. However, the simulated drift rise was sharper than that in our study, which provides a faster drift detection. Baysal et al. [20] detected 88.34% of all real PISA events by using a rule-based method and based on the CGM signal processed by a Kalman filter. Facchinetti et al. [37] modeled PISA as the output of a first-order linear dynamic system driven by a rectangular function, using the data from 72 subjects monitored with

94

Z. Mahmoudi et al. / Biomedical Signal Processing and Control 38 (2017) 86–99

Drift detection flag for sensor 1 Drift of sensor 1

Glucose concentration (mg/dL)

1000 800

800

600

600

400

400

200 100 0 -100 -200

200 100 0 -100 -200

0

4

8 12 16 Time (hour)

20

Drift detection flag for sensor 2 Drift of sensor 2

1000

24

0

4

8 12 16 Time (hour)

20

24

(a)

CGM data of sensor 1 with drift CGM data of sensor 2 with drift Filtered CGM with adaptive UKF Filtered CGM with nonadaptive UKF Simulated BG

Glucose concentration (mg/dL)

1000 800 600 400 200 100 0 -100

0

4

8

12 Time (hour)

16

20

24

20

24

(b) 3

× 10

5

Covariance of the measurement noise for sensor 1 Covariance of the measurement noise for sensor 2

Covariance (mg

2

2

/dL )

2.5 2 1.5 1 0.5 0

0

4

8

12 Time (hour)

16

(c) Fig. 9. Drift detection. (a) Left panel: example of a positive drift in sensor 1, and the detection flag. (a) Right panel: example of a negative drift in sensor 2 and the detection flag. (b) Data of the two CGM sensors with the simulated drifts in (a), accompanied by the simulated BG concentration and the CGM data filtered by the adaptive UKF. For comparison, the CGM data filtered with a nonadaptive UKF is also illustrated. (c) The covariances of the measurement noise for the two sensors, tuned according to the adaptive law in (22d).

the Dexcom G4P CGM. The modeled PISA events had the median of 45 min for the duration and 24 mg/dL for the amplitude. Emami et al. [38] modeled PISA with a second order polynomial. In order to estimate the nadir and the duration of PISA as the parameters of the model, they used real data from 15 type 1 ® diabetes patients using the Sof-sensor glucose sensor (Medtronic, Northridge, CA). Furthermore, they augmented the obtained PISA model with the sensor model of Lunn et al. [39] in order to obtain an enhanced sensor model for the integration with a dual-hormone

AP simulator. The results indicate that the sensor model augmented with the PISA model reduced the deviance information criteria compared to the original Lunn model. The simulated PISA events in our study indicate a reasonable resemblance to the observed real PISA events in the sensor data that we used. The difference between the PISA duration and amplitude of our data and those of the other data in studies such as [20,37] could be due to the use of different CGM sensors. Sensors of different types may react differently to the compression arti-

Z. Mahmoudi et al. / Biomedical Signal Processing and Control 38 (2017) 86–99

95

Glucose concentration (mg/dL)

PISA detection flag for sensor 1 PISA of sensor 1

200

100

0

-100

-200

0

4

8

12 Time (hour)

16

20

24

(a)

Glucose concentration (mg/dL)

PISA detection flag for sensor 2 PISA of sensor 2

200

100

0

-100

-200

0

4

8

12 Time (hour)

16

20

24

(b)

Glucose concentration (mg/dL)

400 CGM data of sensor 1 with PISA CGM data of sensor 2 with PISA Filtered CGM with adaptive UKF Filtered CGM with non-adaptive UKF Simulated BG

300

200

100

0

-100

0

4

8

12 Time (hour)

16

20

24

(c)

Covariance of the measurement noise for sensor 1 Covariance of the measurement noise for sensor 2

8000

Covariance (mg

2

/dlL 2 )

10000

6000

4000

2000

0

0

4

8

12 Time (hour)

16

20

24

(d) Fig. 10. PISA detection. (a) Example of 5 simulated PISA events in sensor 1 and the detection flag. (b) Example of 5 simulated PISA events in sensor 2 and the detection flag. (c) Data of the two CGM sensors with the simulated PISA events in (a) and (b), accompanied by the simulated BG concentration and the CGM data filtered by the adaptive UKF. For comparison, the CGM data filtered with a nonadaptive UKF is also illustrated. (c) The covariance of the measurement noise for the two sensors, tuned according to the adaptive law in (22d).

96

Z. Mahmoudi et al. / Biomedical Signal Processing and Control 38 (2017) 86–99

Glucose concentration (mg/dL)

300 CGM data of sensor 1 CGM data of sensor 2 Meal detection flag Meal taken by the virtual patient Meal announced to the model

200

100

0

0

4

8

12 Time (hour)

16

20

24

16

20

24

(a)

Glucose concentration (mg/dL)

300 CGM data of sensor 1 CGM data of sensor 2 Meal taken by the virtual patient Meal announced to the model Meal detection flag

200

100

0

0

4

8

12 Time (hour)

(b) Fig. 11. An example of meal detection. (a) Unannounced meals. (b) Correctly announced meals.

fact depending on their electronics. Therefore, our method for PISA detection should be tuned to the sensor type when applied to the real data. We simulated drift as a decremental or incremental linear trend between the two calibration instances. Modeling drift as a linear trend appears to fit the real data well, as demonstrated by Emami et al. [38] and Facchinetti et al. [27]. Dassau et al. [40] developed a method for detection of unannounced meals for an AP. Their methods included applying threshold on the glucose rate-of-change by backward difference (BD), the Kalman filter estimation, combination of BD and the Kalman filter, and the second derivative of glucose. To minimize false positive detection they applied a voting algorithm consisting of a two-out-of-three or three-out-of four schemes from the proposed methods. When the method was applied on the data of a FreeStyle Navigator (Abbott Diabetes Care, Alameda, CA), the mean detection delay from the onset of the meal was 30 min and the mean increase in the serum glucose was 15 mg/dL. Lee et al. [41] also developed a method for meal detection and meal size estimation using a Kalman filter approach for an AP with the model predictive control algorithm. Turksoy et al. [42] developed a meal detection algorithm based on a CGM, using the Bergman’s minimal model in an UKF. The method works by comparing the estimated rate of change of appearance of glucose with a threshold. Out of 63 meals for nine real subjects, 61 meals were detected. The average change in glucose levels between the meals and the detection points was 16 mg/dL. The advantage of our proposed meal detection method is that it can not only detect the meal but also distinguish it from the sensor artifact. Other signal processing approaches for CGM fault detection include methods based on wavelet decompo-

sition [43,44], support vector machine [45,46], statistical modeling using historical CGM data [47], rate-limiting filtering of the signal [13,48], and the Kalman filtering [15,49,50]. In this paper, we combine a nonlinear model with an UKF to form predictions for the interstitial glucose concentration measured by two CGMs. Using CGM measurements, we compute the prediction errors and announce the measurement as a faulty measurement if the prediction errors does not comply with a multivariate normal distribution. We test this by checking whether the Hotelling’s test statistics of the prediction errors is outside of the corresponding control region. Zhao et al. [51] addressed the model-based CGM fault detection by means of the principal component analysis (PCA). They built a PCA model of the fault-free sensor glucose data. A new data sample is announced as a fault if the squared prediction error for the residual subspace of the model is outside of a control region. Song et al. [52] used the same PCA-based fault detection in a classification framework. The fusion of the posterior probabilities of the classifier with the output of the PCA fault detectors provided a more accurate fault detection. A comparison between our paper and [51,52] indicates that all three studies use the misfit between a model and the CGM data in order to detect a fault. The models in the three papers describe the fault-free data. The difference between the papers is that we use a nonlinear physiological model describing the physiology of the fault-free data, while Zhao et al. [51] and Song et al. [52] use a linear data-driven PCA model that describes the autocorrelation of the fault-free data. The physiological model in (A.1) has 10 parameters to be estimated for having a subject-specific fault detection. The relatively large number of parameters to be identified can limit the practicality of the method. Alternatively, input-output models with more minimal structures should be considered. An example of such models is the linear model introduced by Boiroux et al. [53,54] that has only two parameters to identify, i.e., the insulin sensitivity and the insulin action time. These parameters are often easily identified and used by patients and clinicians for daily insulin adjustments. Another alternative to the physiological models is the black-box data-driven ARMAX models. Oruklu et al. [55] used an adaptive ARMA model. The model can track the fluctuations in glucose metabolism by applying a change detection method on the recursively identified model parameters. Also, Oruklu et al. [56] suggested using an ARMAX model that can include multiple inputs such as food intake, physical activity, and emotional stimuli. Such a model is used to supplement and enhance the CGM information for the prediction of hyperglycemia and hypoglycemia. A comprehensive review of glucose prediction models is provided in [57]. One of the limitations of our method is that it requires information of meals and insulin infusion which may not be always available with present device technology, e.g., when standalone CGM systems are used. Future device technology may allow the integration of smart insulin pens with CGMs. Such smart pens would be able to store the information of insulin administration and communicate wirelessly with the CGM sensors. For these future integrated systems, the availability of insulin information for CGM fault detection is not out of reach. Our fault detection methodology requires the wear of two sensors and benefits from the assumption of independence between the two CGM sensors. This could be a limitation when the two sensors have the same glucose sensing principles. In order to provide independence in this case, the two sensors should be mounted at two distant locations. This could be practically challenging. With the advances in sensor technology, it would be possible to have two independent sensors in one device, mounted close to each other, provided that the sensors have different sensing principles, e.g., electrochemical and optical sensing.

Z. Mahmoudi et al. / Biomedical Signal Processing and Control 38 (2017) 86–99

Regulating blood glucose in the presence of an unannounced meal is a big challenge for the AP. Therefore, detecting an unannounced meal is crucial for a CGM fault detection module. Due to the principle that we use for the fault isolation, the algorithm can also detect other physiological events such as stress-induced hyperglycemia or carb counting errors (causing overestimation or underestimation of meals). Nevertheless, a physiological event detector, detecting changes in the metabolism from CGM measurements, requires a more sophisticated and multi-layered set of rules in the fault isolator. Another limitation of the current study is that the model used for simulating the patient is the same as the model used in the fault detectors. This situation describes a perfect scenario, while in reality, the models in the UKFs can not have the same level of complexity as that of a patient’s metabolism. In order to make a more realistic simulation, a comprehensive metabolic model should simulate the patient, while the UKFs in the fault detectors should use a minimal model. We used the covariance matching technique to estimate the covariance of the measurement noise. This method has low computational expense and is flexible to accommodate the nonlinear applications [58]. However, it is a suboptimal estimation approach. This method should be compared with the optimization-based noise covariance estimation techniques such as autocovariance least-squares methods [59,60], and maximum likelihood estimation [61,62].

6. Conclusions We presented the basis of a signal processing module for fault detection and fault correction in CGM sensors. In addition, the module is able to detect unannounced meals. We tested the proposed methods in a perfect set up, used as a benchmark. The perfect set up, in terms of models, is when the model in the detectors is the same as the model in the patient simulator which the detectors are tested with. The module consists of two CGM sensors, two fault detectors, a fault isolator, and an adaptive UKF with tunable measurement noise covariance. In comparison to a nonadaptive UKF, the adaptive UKF reduces the deviation of the CGM measurements from their paired BG concentrations. Both features, i.e., fault detection/correction and meal detection, are fundamental for an AP application. Fault detection improves the CGM reliability and meal detection paves the way for the realization of a fully automated AP. Another possible structure for the adaptive UKF is tuning the process noise covariance, Q, when the fault isolator detects an unannounced meal or an unforeseen change in the patient’s dynamics. Using this adaptive scheme in the AP allows an adjustable uncertainty in the process model which provides a more robust closed-loop control of the BG concentrations. Further studies are required to evaluate the method using real CGM data.

Funding source This work is funded by the Danish Diabetes Academy supported by the Novo Nordisk Foundation. The funding sources were not involved in the preparation of the manuscript.

97

Appendix A. The virtual patient model The equation set (A.1) describes a SDE model of a type 1 diabetes patient [23]. We use this model in simulation, and in the UKFs for prediction and filtering. 1 1

dI SC (t) =

 u(t) CI



− ISC (t) dt + SC · dωSC (t),

(A.1a)

 1  ISC (t) − Ip (t) dt + p · dωp (t), 2

dI p (t) =

dI eff (t) =





−P2 · Ieff (t) + P2 · SI · Ip (t) dt

+eff · dωeff (t),



dGB (t) =

(A.1b)

(A.1c)



−(GEZI + Ieff (t)) · GB (t) + EGP + RA (t) dt

(A.1d)

+G · dωG (t), dGI (t) = − dd1 (t) = dd2 (t) = RA (t) =



1 (GB (t) − GI (t)) dt + GI · dωGI (t), 3

(A.1e)



1 d1 (t) dt + d1 · dωd1 (t), m

(A.1f)

1 (d1 (t) − d2 (t)) dt + d2 · dωd2 (t), m

(A.1g)

d(t) −

1 d2 (t). VG m

(A.1h)

The parameter  1 (49 min) is the time constant of the insulin movement from the administration site to the SC tissue,  2 (47 min) is the time constant related to the insulin movement from the SC tissue to plasma,  3 (10 min) is the time constant of the glucose movement from plasma to SC tissue, CI (2010 mL/min) is the insulin clearance rate, P2 (1.06 · 10−2 min−1 ) is the delayed insulin action on the GB concentration, SI (8.11 · 10−4 mL/␮U/min) is the insulin sensitivity, GEZI (2.20 · 10−3 min−1 ) is the glucose effectiveness at zero insulin, EGP (1.33 mg/dL/min) is the endogenous glucose production rate at zero insulin,  m (47 min) is the peak time of meal absorption, and Vg (253 dL) is the volume of distribution for glucose. We simulated the model in (A.1) using Euler Maruyama method [63]. The patient eats breakfast at 8:00 h, lunch at 13:15 h, dinner at 18:00 h, and snack at 22:00 h. The CHO content of the meals are 72 g for breakfast, 131 g for lunch, 51 g for dinner, and 70 g for snack, respectively. The insulin boluses are 3.5 U for breakfast, 7 U for lunch, 2.5 U for dinner, and 3.5 U for snack, respectively. Appendix B. Filter parameters (i)

In the fault detectors (local UKFs), Qn =  T , and R = Rn for i = 1, (1) (2) 2. We set Rn = Rn = Rn , and  = 0.5% xss , where xss is the steady (1) (2) state of the model in Appendix A. In general, Rn and Rn can be different. T  In the main UKF, Qn =  , where  = 0.5% xss , and R = R(1) 0 . R(2) 0 For simulating the CGM data for a patient with the model in (A.1),  = 0.5% xss is used. Then n in (2) is added to the simulated GI , and depending on which fault to be detected, drift or PISA is added to the CGM data. References

Conflict of interest The authors declare that there are no conflicts of interest.

[1] B.W. Bequette, Fault detection and safety in closed-loop artificial pancreas systems, J. Diabetes Sci. Technol. 8 (6) (2014) 1204–1214. [2] B.W. Bequette, Challenges and recent progress in the development of a closed-loop artificial pancreas, Annu. Rev. Control 36 (2) (2012) 255–266.

98

Z. Mahmoudi et al. / Biomedical Signal Processing and Control 38 (2017) 86–99

[3] B.W. Bequette, N. Baysal, D. Howsmon, F. Cameron, Steps towards a closed-loop artificial pancreas, 41st Annual Northeast Biomedical Engineering Conference (NEBEC) (2015) 1–2. [4] F. Doyle, L. Huyett, J. Lee, H. Zisser, E. Dassau, Closed-loop artificial pancreas systems: engineering the algorithms, Diabetes Care 37 (5) (2014) 1191–1197. [5] E. Zschornack, C. Schmid, S. Pleus, M. Link, H.-M. Klötzer, K. Obermaier, M. Schoemaker, M. Strasser, G. Frisch, G. Schmelzeisen-Redeker, C. Haug, G. Freckmann, Evaluation of the performance of a novel system for continuous glucose monitoring, J. Diabetes Sci. Technol. 7 (4) (2013) 815–823. [6] T. Bailey, B.W. Bode, M.P. Christiansen, L.J. Klaff, S. Alva, The performance and usability of a factory-calibrated flash glucose monitoring system, Diabetes Technol. Ther. 17 (11) (2015) 787–794. [7] L. Laffel, Improved accuracy of continuous glucose monitoring systems in pediatric patients with diabetes mellitus: results from two studies, Diabetes Technol. Ther. 18 (2) (2016) 23–33. [8] J. Kropff, J.H. DeVries, Continuous glucose monitoring, future products, and update on worldwide artificial pancreas projects, Diabetes Technol. Ther. 18 (2) (2016) 53–63. [9] A. Facchinetti, Continuous glucose monitoring sensors: past, present and future algorithmic challenges, Sensors 16 (12) (2016) 2093. [10] J.R. Castle, P.G. Jacobs, Nonadjunctive use of continuous glucose monitoring for diabetes treatment decisions, J. Diabetes Sci. Technol. 10 (5) (2016) 1169–1173. [11] W.K. Ward, A review of the foreign-body response to subcutaneously-implanted devices: the role of macrophages and cytokines in biofouling and fibrosis, J. Diabetes Sci. Technol. 2 (5) (2008) 768–777. [12] K.L. Helton, B.D. Ratner, N.A. Wisniewski, Biomechanics of the sensor-tissue interface-effects of motion, pressure, and design on sensor performance and foreign body response – part II: examples and application, J. Diabetes Sci. Technol. 5 (3) (2011) 647–656. [13] Z. Mahmoudi, M. Dencker-Johansen, J. Christiansen, O. Hejlesen, A multistep algorithm for processing and calibration of microdialysis continuous glucose monitoring data, Diabetes Technol. Ther. 15 (10) (2013) 825–835. [14] G.M. Steil, K. Rebrin, P.V. Goode, jr., J.J. Mastrototaro, R.R. Purvis, W.P. Van Antwerp, J.J. Shin, C.D. Talbot, Closed loop system for controlling insulin infusion. U.S. patent number: CA 2373986 C (2005). [15] A. Facchinetti, G. Sparacino, C. Cobelli, Online denoising method to handle intraindividual variability of signal-to-noise ratio in continuous glucose monitoring, IEEE Trans. Biomed. Eng. 58 (8) (2011) 2664–2671. [16] A. Facchinetti, S. Del Favero, G. Sparacino, C. Cobelli, Detecting failures of the glucose sensor-insulin pump system: improved overnight safety monitoring for type-1 diabetes, in: Annual International Conference of the IEEE Engineering in Medicine and Biology Society, IEEE, 2011, pp. 4947–4950. [17] J. Feng, K. Turksoy, S. Samadi, I. Hajizadeh, A. Cinar, Hybrid online sensor error detection and functional redundancy for artificial pancreas control systems, Proceeding of the 11th IFAC Symposium on Dynamics and Control of Process Systems, including Biosystems (2016) 753–758. [18] H. Zhao, C. Zhao, An automatic denoising method with estimation of noise level and detection of noise variability in continuous glucose monitoring, Proceeding of the 11th IFAC Symposium on Dynamics and Control of Process Systems, including Biosystems (2016) 785–790. [19] K. Turksoy, L. Quinn, E. Littlejohn, A. Cinar, Monitoring and fault detection of continuous glucose sensor measurements, Proceedings of the American Control Conference (2015) 5091–5096. [20] N. Baysal, F. Cameron, B.A. Buckingham, D.M. Wilson, H.P. Chase, D.M. Maahs, B.W. Bequette, A novel method to detect pressure-induced sensor attenuations (PISA) in an artificial pancreas, J. Diabetes Sci. Technol. 8 (6) (2014) 1091–1096. [21] Z. Mahmoudi, D. Boiroux, M. Hagdrup, K. Nørgaard, N.K. Poulsen, H. Madsen, J.B. Jørgensen, Application of the continuous-discrete extended Kalman filter for fault detection in continuous glucose monitors for people with type 1 diabetes, Proceedings of the European Control Conference (2016) 714–719. [22] Z. Mahmoudi, S.L. Wendt, D. Boiroux, M. Hagdrup, K. Nørgaard, N.K. Poulsen, H. Madsen, J.B. Jørgensen, Comparison of three nonlinear filters for fault detection in continuous glucose monitors, Proceedings of the 38th Annual International Conference of the IEEE Engineering in Medicine and Biology Society (2016) 3507–3510. [23] S. Kanderian, S. Weinzimer, G. Voskanyan, G. Steil, Identification of intraday metabolic profiles during closed-loop glucose control in individuals with type 1 diabetes, J. Diabetes Sci. Technol. (2009) 1047–1057. [24] R. Hovorka, L.J. Chassin, M. Ellmerer, J. Plank, M.E. Wilinska, A simulation model of glucose regulation in the critically ill, Physiol. Meas. 29 (8) (2008) 959–978. [25] M.E. Wilinska, L.J. Chassin, C.L. Acerini, J.M. Allen, D.B. Dunger, R. Hovorka, Simulation environment to evaluate closed-loop insulin delivery systems in type 1 diabetes, J. Diabetes Sci. Technol. 4 (1) (2010) 132–144. [26] L.J. Chassin, M.E. Wilinska, R. Hovorka, Evaluation of glucose controllers in virtual environment: methodology and sample application, Artif. Intell. Med. 32 (3) (2004) 171–181. [27] A. Facchinetti, S. Del Favero, G. Sparacino, J. Castle, W. Ward, C. Cobelli, Modeling the glucose sensor error, IEEE Trans. Biomed. Eng. 61 (3) (2014) 620–629. [28] D. Simon, Optimal State Estimation: Kalman, H∞, and Nonlinear Approaches, John Wiley & Sons, 2006.

[29] S. Särkkä, On unscented Kalman filtering for state estimation of continuous-time nonlinear systems, IEEE Trans. Autom. Control 52 (9) (2007) 1631–1641. [30] S. Julier, J. Uhlmann, Unscented filtering and nonlinear estimation, Proc. IEEE 92 (3) (2004) 401–422. [31] J.B. Jørgensen, S.B. Jørgensen, Comparison of prediction-error-modelling criteria, Proceedings of the American Control Conference (2007) 5300–5306. [32] S.X. Ding, Data-driven Design of Fault Diagnosis and Fault-tolerant Control Systems, Springer London, 2014. [33] G. Li, S.J. Qin, Comparative study on monitoring schemes for non-Gaussian distributed processes, J. Process Control (2016), http://dx.doi.org/10.1016/j. jprocont.2016.08.007. [34] S.J. Qin, Survey on data-driven industrial process monitoring and diagnosis, Annu. Rev. Control 36 (2) (2012) 220–234. [35] P.S. Maybeck, Stochastic Models, Estimation, and Control, Academic Press, 1982 (Chapter 10). [36] A. Savitzky, M.J. Golay, Smoothing and differentiation of data by simplified least squares procedures, Anal. Chem. 36 (1964) 1627–1639. [37] A. Facchinetti, S. Del Favero, G. Sparacino, C. Cobelli, Modeling transient disconnections and compression artifacts of continuous glucose sensors, Diabetes Technol. Ther. 18 (4) (2016) 264–272. [38] A. Emami, R. Rabasa-Lhoret, A. Haidar, Enhancing glucose sensor models: modeling the drop-outs, Diabetes Technol. Ther. 17 (6) (2015) 420–426. [39] D. Lunn, C. Wei, R. Hovorka, Fitting dynamic models with forcing functions: application to continuous glucose monitoring in insulin therapy, Stat. Med. 30 (18) (2011) 2234–2250. [40] E. Dassau, B.W. Bequette, B.A. Buckingham, F.J. Doyle, Detection of a meal using continuous glucose monitoring implications for an artificial beta-cell, Diabetes Care 31 (2) (2008) 295–300. [41] H. Lee, B.A. Buckingham, D.M. Wilson, B.W. Bequette, A closed-loop artificial pancreas using model predictive control and a sliding meal size estimator, J. Diabetes Sci. Technol. 3 (5) (2009) 1082–1090. [42] K. Turksoy, S. Samadi, J. Feng, E. Littlejohn, L. Quinn, A. Cinar, Meal detection in patients with type 1 diabetes: a new module for the multivariable adaptive artificial pancreas control system, IEEE J. Biomed. Health Inform. 20 (1) (2016) 47–54. [43] Q. Shen, S. Qin, K. Doniger, Online dropout detection in subcutaneously implanted continuous glucose monitoring, Proceedings of the American Control Conference (2010) 4373–4378. [44] G. Campetelli, D. Zumoffen, M. Basualdo, Improvements on noninvasive blood glucose biosensors using wavelets for quick fault detection, J. Sens. (2011) 1–11. [45] Y. Leal, L. Gonzalez-Abril, C. Lorencio, J. Bondia, J. Vehi, Detection of correct and incorrect measurements in real-time continuous glucose monitoring systems by applying a postprocessing support vector machine, IEEE Trans. Biomed. Eng. 60 (7) (2013) 1891–1899. [46] J. Bondia, C. Tarín, W. García-Gabin, E. Esteve, J. Fernández-Real, W. Ricart, J. Vehí, Using support vector machines to detect therapeutically incorrect ® measurements by the MiniMed CGMS , J. Diabetes Sci. Technol. 2 (4) (2008) 622–629. [47] M. Signal, A. Le Compte, D. Harris, P. Weston, J. Harding, J. Chase, Using stochastic modelling to identify unusual continuous glucose monitor measurements and behaviour, in newborn infants, BioMed. Eng. OnLine 11 (1) (2012) 45. [48] N. Baysal, F. Cameron, B.A. Buckingham, D.M. Wilson, B.W. Bequette, Detecting sensor and insulin infusion set anomalies in an artificial pancreas, Proceedings of the American Control Conference (2013) 2929–2933. [49] S. Del Favero, M. Monaro, A. Facchinetti, A. Tagliavini, G. Sparacino, C. Cobelli, Real-time detection of glucose sensor and insulin pump faults in an artificial pancreas, Proceedings of the 19th IFAC World Congress (2014) 1941–1946. [50] A. Facchinetti, S. Del Favero, G. Sparacino, C. Cobelli, An online failure detection method of the glucose sensor-insulin pump system: improved overnight safety of type-1 diabetic subjects, IEEE Trans. Biomed. Eng. 60 (2) (2013) 406–416. [51] C. Zhao, Y. Fu, Statistical analysis based online sensor failure detection for continuous glucose monitoring in type I diabetes, Chemom. Intell. Lab. Syst. 144 (2015) 128–137. [52] G. Song, C. Zhao, Y. Sun, A classification-based fault detection method for continuous glucose monitoring (CGM), Proceedings of the World Congress on Intelligent Control and Automation (WCICA) (2016) 956–961. [53] D. Boiroux, A.K. Duun-Henriksen, S. Schmidt, K. Nørgaard, N. Poulsen, H. Madsen, J.B. Jørgensen, Adaptive control in an artificial pancreas for people with type 1 diabetes, Control Eng. Pract. 58 (20) (2017) 332–342. [54] D. Boiroux, V. Bátora, M. Hagdrup, M. Tárnik, J. Murgaˇs, S. Schmidt, K. Nørgaard, N.K. Poulsen, H. Madsen, J.B. Jørgensen, Comparison of prediction models for a dual-hormone artificial pancreas, IFAC Workshop Ser. 48 (20) (2015) 7–12. [55] M. Eren-Oruklu, A. Cinar, L. Quinn, D. Smith, Adaptive control strategy for regulation of blood glucose levels in patients with type 1 diabetes, J. Process Control 19 (8) (2009) 1333–1346. [56] M. Eren-Oruklu, A. Cinar, D.K. Rollins, L. Quinn, Adaptive system identification for estimating future glucose concentrations and hypoglycemia alarms, Automatica 48 (8) (2012) 1892–1897. [57] Prediction Methods for Blood Glucose Concentration, in: H. Kirchsteiger, J.B. Jørgensen, E. Renard, L. Del Re (Eds.), Springer, 2016.

Z. Mahmoudi et al. / Biomedical Signal Processing and Control 38 (2017) 86–99 [58] K.A. Myers, B.D. Tapley, Adaptive sequential estimation with unknown noise statistics, IEEE Trans. Autom. Control AC-21 (4) (1976) 520–523. [59] B.J. Odelson, A. Lutz, J.B. Rawlings, The autocovariance least-squares method for estimating covariances: application to model-based control of chemical reactors, IEEE Trans. Control Syst. Technol. 14 (3) (2006) 532–540. [60] M.A. Zagrobelny, J.B. Rawlings, Practical improvements to autocovariance least-squares, AIChE 61 (6) (2015) 1840–1855.

99

[61] N.R. Kristensen, H. Madsen, S.B. Jørgensen, Parameter estimation in stochastic grey-box models, Automatica 40 (2) (2004) 225–237. [62] M.A. Zagrobelny, J.B. Rawlings, Identifying the uncertainty structure using maximum likelihood estimation, Proc. Am. Control Conference 2015 (2015) 422–427. [63] D.J. Higham, An algorithmic introduction to numerical simulation of stochastic differential equations, SIAM Rev. 43 (3) (2001) 525–546.