Sage Journals: Discover world-class research

Abstract

Background:

One way of constructing a control algorithm for an artificial pancreas is to identify a model capable of predicting plasma glucose (PG) from interstitial glucose (IG) observations. Stochastic differential equations (SDEs) make it possible to account both for the unknown influence of the continuous glucose monitor (CGM) and for unknown physiological influences. Combined with prior knowledge about the measurement devices, this approach can be used to obtain a robust predictive model.

Method:

A stochastic-differential-equation-based gray box (SDE-GB) model is formulated on the basis of an identifiable physiological model of the glucoregulatory system for type 1 diabetes mellitus (T1DM) patients. A Bayesian method is used to estimate robust parameters from clinical data. The models are then used to predict PG from IG observations from 2 separate study occasions on the same patient.

Results:

First, all statistically significant diffusion terms of the model are identified using likelihood ratio tests, yielding inclusion of $σ_{I_{s c}}$ , $σ_{G_{p}}$ , and $σ_{G_{s c}}$ . Second, estimates using maximum likelihood are obtained, but prediction capability is poor. Finally a Bayesian method is implemented. Using this method the identified models are able to predict PG using only IG observations. These predictions are assessed visually. We are also able to validate these estimates on a separate data set from the same patient.

Conclusions:

This study shows that SDE-GBs and a Bayesian method can be used to identify a reliable model for prediction of PG using IG observations obtained with a CGM. The model could eventually be used in an artificial pancreas.

Keywords

Bayesian methods plasma glucose dynamics PG-IG dynamics stochastic differential equations stochastic gray-box modeling type 1 diabetes mellitus

The use of continuous glucose monitors (CGMs) on patients with type 1 diabetes mellitus (T1DM) has a positive effect on glycemic control.^1-5 The advantages of continuous glucose monitoring are many, but there are also several well-documented problems.^6-8 These include time delay between plasma glucose (PG) and interstitial glucose (IG),^3-5 the drift of sensor sensitivity,⁹ and vulnerability to calibrations during rapidly increasing or decreasing glucose levels.^10,11 For optimal use of CGMs in the management of T1DM, for example, for closed-loop control, it is necessary to know the relationship between PG and IG.¹² Due to the well-known issues with the sensor, the main challenge related to automatic control of the glucose is to predict PG from IG observations measured by a CGM. To do this, the dynamical relation between the PG and IG levels needs to be well described.

Modeling glucose dynamics using observations from CGMs is challenging. As described by Facchinetti et al., three of the major problems are (1) PG observations have to be collected parallel with IG observations, (2) a precise model for PG-IG dynamics is needed, and (3) sensor errors should be modeled as stochastic processes.¹³

When using real-life data obtained with a CGM, the autocorrelation of data becomes an issue, requiring specific modeling techniques. This is discussed in this article. As the observations are closely sampled (e.g., 5-minute sample rate), interdependence between data points or autocorrelation is present.^14-16 Several approaches to overcome sensor issues, including autocorrelation, have been investigated.^17-19 Data-driven autoregressive models have been used by Reifman et al.¹⁶ Chase et al. modeled sensor errors using Gaussian processes.²⁰ Breton and Kovatchev used a non-Gaussian Johnson distribution and autoregressive (AR(1) processes).²¹ One suggestion by Guerra et al. is to use deconvolution-based enhancement algorithms that, in real time, improve CGM accuracy.²²

In this study we use stochastic-differential-equation-based gray-box (SDE-GB) models and a Bayesian method to address all 3 issues described by Facchinetti et al.

When investigating the glucose-insulin dynamics (e.g., in closed-loop studies), patients receive inputs to excitate glucose values, such as meals, giving a positive excitation and insulin, giving a negative excitation. This article accommodates the need for an understanding of the dynamics with no external positive excitation. The negative excitation is provided by a continuous subcutaneous insulin infusion (CSII).

Glucose-insulin dynamics are a part of the complex, partly unknown fuel metabolism.²³ Stochastic behavior should therefore be considered in the modeling. One way of doing this is to apply SDE-GB models. Such models are established using a combination of physiological knowledge and statistical information from data. In general a SDE-GB model can be written as,

d x_{t} = f (x_{t}, u_{t}, t, θ) d t + σ (x_{t}, u_{t}, t, θ) d ω

(1)

y_{j} = h (x_{j}, u_{k}, t_{j}, θ) + e_{j}

(2)

where the states are described by the vector $x_{t}$ , $f (\cdot)$ is the drift term and $σ (\cdot) d ω$ is a diffusion term with $ω$ as a Wiener process.²⁴ The states are unobserved, but the vector, $y_{j},$ of discrete time observations is linked to the states as described by the function $h (\cdot) .$ The observations are contaminated by the measurement noise, $e_{j} .$ The input is represented by $u_{t}$ and $θ$ represents the parameters.

The advantage of using SDE-GB models as compared to ordinary differential-equations-based models, is that noise is included in the state equations in terms of $σ (\cdot) d ω .$ As the noise affects the state equations, the overall system will have a stochastic influence, reflecting the fact that future values of the states can never be predicted exactly. The system noise enters through the dynamics of the system and will reveal, for example, model deficiencies, random fluctuations in the system or unknown inputs.²⁵ Residual information and analysis can be used in the validation process of the model building to ensure that uncorrelated residuals are obtained.²⁴ Uncorrelated residuals indicate that the model is adequate in the sense that all the systematic variation in the observations is described. Elaborations of these advantages exist in several articles.^25-30 Solutions to SDEs are stochastic processes. These stochastic processes are described by probability distributions.³¹ This property enables the use of maximum likelihood techniques, thereby maintaining a data-driven model estimation.³² In this article, the parameters of the Medtronic virtual patient (MVP)³³ model are estimated in an SDE-GB setting using the maximum likelihood principle. First, a method for identifying the need for diffusion terms, using forward selection, is presented. Second, parameters are estimated using maximum likelihood. Finally a Bayesian method is used to estimate the parameters related to the observation equations. These parameters are cross-validated on a different data set from the same patient.

Methods

The data used in this study originate from an overnight closed-loop study, conducted as a part of the DIACON project.³⁴ The patient was studied on 2 overnight hospital admissions from 17:30 to 07:00. At 18:00, an evening meal was served. The size of the meal was determined by patient weight (1 g carbohydrate/kg body weight) and was identical on the 2 study nights. The bolus calculator in the patient’s own insulin pump was used to determine the size of the meal bolus. At 22:00 the study initiated. No further food or insulin boluses were given.³⁴ PG and plasma insulin were sampled every 30 minutes. PG was analyzed using YSI2300 STAT Plus (Yellow Springs Instruments, Yellow Springs, OH). Every 5 minutes a subcutaneous glucose value was obtained using a Dexcom SEVEN PLUS (Dexcom, San Diego, CA) CGM. Insulin Aspart (Novorapid, Novo Nordisk, Bagsværd, Denmark) was administered using a Paradigm Veo (Medtronic, Northridge, CA) insulin pump. Patients were euglycemic (PG level 70-144 mg/dl) at study start. It is assumed that an evening meal given 4 hours before study start does not influence the glucose dynamics. For this study, data from on 1 patient on 2 independent nights at least 2 weeks apart were used. The first data set, Data₁, originates from a closed-loop study night, and the second data set, Data₂, originates from an open-loop study night.

The Initial Model

The MVP model was used as the basis for formulating the SDE-GB model. In essence, the model is based on the Bergman minimal model.^35,36 Fisher modified the model, replacing insulin secretion with insulin infusion.³⁷ In this study, we focus on the kinetics and dynamics when patients are at rest. This reduces the model to the 5 equations below describing glucose-insulin dynamics with no meal input.

The insulin dynamics are described by,

d I_{s c} = \frac{1}{τ_{1}} (\frac{I D}{C} - I_{s c}) d t + σ_{I_{s c}} d ω_{1}

(3)

d I_{P} = \frac{1}{τ_{2}} (I_{s c} - I_{P}) d t + σ_{I_{P}} d ω_{2}

(4)

d I_{e f f} = (p_{2} \cdot S_{I} \cdot I_{P} - p_{2} \cdot I_{e f f}) d t + σ_{I_{e f f}} d ω_{3}

(5)

where $I_{s c}, I_{P}$ and $I_{e f f}$ are the subcutaneous insulin concentration $[\frac{I U}{L}]$ , the plasma insulin concentration $[\frac{I U}{L}]$ and the effect of insulin $[\frac{1}{m i n}],$ respectively.

The glucose dynamics are described as,

d G_{P} = (- (G E Z I + I_{e f f}) \cdot G_{P} + E G P) d t + σ_{G_{P}} d ω_{4}

(6)

d G_{s c} = \frac{1}{τ_{3}} (- G_{s c} + G_{P}) d t + σ_{G_{s c}} d ω_{5}

(7)

where $G_{P}$ and $G_{s c}$ are states representing the PG concentration $[\frac{mg}{dl}]$ and subcutaneous IG concentration $[\frac{mg}{dl}]$ , respectively. The models contain a diffusion term $σ d ω$ where $σ$ is a scaling parameter and $ω$ is assumed to be a Wiener process with normally distributed increments.²⁸ To accommodate available data, in our study $τ_{2} = τ_{1},$ as insulin kinetics is very data dependent.³⁸ The remaining parameters are defined in Table 1. The model can be expanded to include extra observations equations and state equations accounting for, for example, fingerstick calibrations.

Table 1.

Identifiable Parameters in Medtronic Virtual Patient on SDE-GB Form.

Name	Unit	Description
EGP	mg/(dL×min)	Endogenous glucose production rate, estimated at zero insulin
GEZI	1/min	Effect of glucose per se to increase glucose uptake into cells and lower endogenous glucose production at zero insulin
$p_{2}$	1/min	The delay in insulin action following an increase in plasma insulin
$S_{i}$	L/(mU×min)	Insulin sensitivity
$σ_{I_{S C}}$	—	Scaling of diffusion term for interstitial insulin
$σ_{I_{P}}$	—	Scaling of diffusion term for blood insulin^a
$σ_{I_{e f f}}$	—	Scaling of diffusion term for insulin effect^a
$σ_{G_{P}}$	—	Scaling of diffusion term for PG observation
$σ_{G_{S C}}$	—	Scaling of diffusion term for IG observation
$S_{Y S I}$	—	Variance of measurement noise of the PG observations
$S_{C G M}$	—	Variance of measurement noise of the IG observations
$S_{I N S}$	—	Variance of measurement noise of the insulin observations
$τ_{1}$	min	Time constant related to the insulin movement between the subcutaneous tissue and plasma
$τ_{3}$	min	Time constant related to the glucose movement in the subcutaneous layer
ID	IU/min	INPUT. Insulin delivery to the subcutaneous layer
C	L/min	Clearance of insulin in the subcutaneous tissue (1.2)^b
BA	-	Bioavailability of insulin (0.7)^c

IG, interstitial glucose; PG, plasma glucose.

Found to be insignificant—see results.

Clearance is fixed conservatively according to patient characteristic.^32-34

Bioavailability is fixed conservatively according to Bergman et al.³⁵

For our study the observation equations are,

Y S I = G_{P} + e x p (e_{Y S I})

(8)

C G M = G_{s c} + e x p (e_{C G M})

(9)

I N S = I_{P} + \exp (e_{I N S})

(10)

where $e_{Y S I} \in N (0, S_{Y S I})$ , $e_{C G M} \in N (0, S_{C G M})$ and $e_{I N S} \in N (0, S_{I N S})$ are mutually independent white noise processes, and S is the variance of each measurement noise, respectively.

The Maximum Likelihood Principle

The likelihood function captures all the information in the data about the parameters. For time-series data, the likelihood function is a product of 1-step conditional densities. Due to the Gaussian increments of the Wiener processes, it is assumed that the conditional densities are Gaussian. To evaluate this, 1-step conditional predictions are needed. These are calculated using an extended Kalman filter,^26,27 and the likelihood function has been identified as the key inferential quantity conveying all information in statistical modeling including uncertainty.^32,39

Using the 1-step prediction error, $ϵ_{k},$ and the associated variances, $R_{k | k - 1},$ for the extended Kalman filter, it is possible to formulate a Bayesian method. Bayes’s theorem can be applied to obtain an improved estimate by forming the posterior probability density function,³²

p (θ | Υ_{N}) = \frac{p (Υ_{N} | θ) p (θ)}{p (Υ_{N})} \propto p (Υ_{N} | θ) p (θ)

(11)

where $Υ_{N}$ is the set of observations and $θ$ is a vector of parameters.

Assuming Gaussian distributions, the posterior probability density function becomes,

\begin{array}{l} p (θ | Υ_{N}) = (\prod_{k = 1}^{N} \frac{e x p (\frac{1}{2} ϵ_{k}^{T} R_{k | k - 1}^{- 1} ϵ_{k})}{\sqrt{\det (R_{k | k - 1})} {(\sqrt{2 π})}^{\dim (Υ_{N})}}) p (y_{0} | θ) \\ \times p r i o r (θ_{P}, Σ_{P}) \end{array}

(12)

where $y_{0}$ is the initial conditions and $θ_{P}$ and $Σ_{P}$ are the mean and covariance of a prior distribution. The parameter estimates are found by minimizing the negative log-posterior probability density function,

\hat{θ} = \arg \min_{} {- l n (p (θ | Υ_{N}, y_{0}))}

(13)

where $\hat{θ}$ is the Bayesian estimate. This estimate and the covariance matrix, can be used as new prior knowledge for further estimation.

To mimic potential future applications, 3 different setups for the prior knowledge have been considered:

A theoretical prior knowledge for the variance of the measuring devices: $p r i o r (θ_{M D})$ , shown in the appendix. The prior knowledge of the YSI device is based on Khan et al.⁴⁰ For the Dexcom SEVEN PLUS it is set as arbitrary but based on the user guide.⁴¹ For the insulin measurement device it is set as arbitrary.

A prior knowledge formed by the estimate and covariance $(θ_{11}, Σ_{θ_{11}})$ of Data₁: $p r i o r (θ_{11})$ .

A prior knowledge formed by the estimate and covariance $(θ_{22}, Σ_{θ_{22}})$ of Data₂: $p r i o r (θ_{22})$ .

Modeling Process

All computations are performed using the statistical software R (version 2.15.1) and the CTSM-R package.⁴² The identified models are cross-validated by analyzing the predictions from 2 study nights from the same patient. The following steps are followed for both data sets (see Table 2).

Table 2.

Flow of the Modeling Process.

	Step 1: ML estimation^a	Step 2: Bayesian estimation^a using $p r i o r (θ_{M D})$	Step 3: Prediction^b	Step 4: Form prior	Step 5: Cross-validation Bayesian estimation^a	Step 6: Prediction^b
Data₁	${\hat{θ}}_{1}$	${\hat{θ}}_{11}$	${\hat{X}}_{1}$	Prior(θ₁₁)	${\hat{θ}}_{21}$	${\hat{X}}_{21}$
Data₂	${\hat{θ}}_{2}$	${\hat{θ}}_{12}$	${\hat{X}}_{2}$	Prior(θ₁₂)	${\hat{θ}}_{22}$	${\hat{X}}_{22}$

IG, interstitial glucose; PG, plasma glucose.

All observations in the data set are used for estimation.

Only IG observations are used in the prediction of PG in this step.

Step 1 (Diffusion Term Identification and Maximum Likelihood)

In the initial model building process for diffusion step inclusion, forward selection and Likelihood ratio tests are used.^26,27 (a) The parameters of the initial ordinary differential equation model, $σ = 0$ , are estimated. (b) The diffusion terms are included, 1 at a time, and the parameters are estimated again. (c) Wilks’s likelihood ratio test is applied to identify the most significant improvement.³² (d) Steps b-c are repeated until addition of extra diffusion terms does not give a statistically significant improvement. (e) For Data₁ and Data₂ each respective maximum likelihood parameter set $({\hat{θ}}_{1}, Σ_{{\hat{θ}}_{1}})$ and $({\hat{θ}}_{2}, Σ_{{\hat{θ}}_{2}})$ is estimated.

Steps 2-6 (A Bayesian Method)

Step 2

For Data₁ and Data₂, each respective Bayesian estimate, $({\hat{θ}}_{11}, Σ_{{\hat{θ}}_{11}})$ and $({\hat{θ}}_{22}, Σ_{{\hat{θ}}_{22}})$ , is obtained using the measurement device prior knowledge $p r i o r (θ_{M D})$ .

Step 3

With the estimates obtained in step 2, the predictions of PG for Data₁ and Data₂, are calculated using only the respective IG to verify the prediction capabilities.

Step 4

The new prior knowledge $p r i o r (θ_{11})$ and $p r i o r (θ_{22})$ is formed from $({\hat{θ}}_{11}, Σ_{{\hat{θ}}_{11}})$ and $({\hat{θ}}_{22}, Σ_{{\hat{θ}}_{22}})$ .

Step 5

From Data₁, $p r i o r (θ_{11})$ , is used in a new Bayesian estimation using Data₂. Thereby a cross-validation estimation, giving a new Bayesian parameter estimation ${\hat{θ}}_{12}$ (subscript is read as: prior knowledge from Data₁ used to estimate Data₂) is obtained. Likewise with $p r i o r (θ_{22})$ .

Step 6

With the obtained estimates ${\hat{θ}}_{21}$ and ${\hat{θ}}_{12}$ , the predictions of PG for Data₁ and Data₂ are obtained using the respective IG observations only.

Results

For all estimates presented, identifiability has been established by testing different starting points of the estimation. The forward selection method gave statistical reason for diffusion term inclusion and resulted in the terms $σ_{I_{s c}}$ , $σ_{G_{p}}$ and $σ_{G_{s c}}$ . For all results involving Data₁ an outlier was removed in the estimation and prediction. It is represented by a purple circle in the top plot in the relevant figures. The first PG observation is used as a calibration point.

Maximum Likelihood (Step 1)

For the maximum likelihood investigation only the results from Data₁ are presented. Results from Data₂ show identical modeling behavior. Using CTSM-R, the 1-step-ahead prediction of PG and IG is obtained with calibration at each data point. As seen in Figure 1, the maximum likelihood estimate, $θ_{1}$ , will seek to explain all the variation in the data without reflecting that the truth is better described by the PG observations. Thus the prediction of PG shows the unwanted IG dynamics.

Figure 1.

(Overall) Predicting with ${\hat{θ}}_{1}$ using both PG and IG observations. The green ribbon reflects the normoglycemic range. The purple observation in the top plot is an outlier not used in the maximum likelihood estimation nor prediction. (Top) One-step prediction and 95% prediction band of the PG observations. The prediction is not in agreement with the observations. Steep transitions are seen when there is a change in the IG observations. (Middle) One-step prediction of the IG observations and 95% prediction band. The prediction band is unrealistically narrow, reflecting the fact that the model seeks to explain the sensor error exactly, not reflecting that the PG observations are the truth. (Bottom) Insulin delivery as boluses.

It is evident that the overall prediction of PG follows the dynamic of the IG observations. Steep transitions in the prediction of PG are a direct result of a steep transition in IG observations. These transitions and the overall dynamic behavior are attributed to the maximum likelihood estimate, visualizing the challenge of the estimation technique. The parameter estimates are not presented, as the visual inspection shows suboptimal predictions. This is most easily seen by the persistent deviation between the predictions and the data over longer periods for the PG observations.

A Bayesian Method (Step 2)

Figure 2 shows the predictions from a Bayesian estimation, ${\hat{θ}}_{11}$ , by implementation of prior knowledge, $p r i o r (θ_{M D})$ . It is evident that the prediction of PG follows the actual PG observations more precisely than Figure 1.

Figure 2.

(Overall) Prediction with ${\hat{θ}}_{11}$ using both PG and IG observations. Green ribbon reflects the normoglycemic range. The purple observation in the top plot is an outlier not used in the Bayesian estimation nor prediction. (Top) One-step prediction and 95% prediction band of the PG observations. Small adjustments are seen when a new PG observation is obtained, but the prediction is in agreement with the observations. (Bottom) One-step prediction of IG observations and 95% prediction band.

The prediction of IG is not as precise as before and the confidence interval is wider, a direct consequence of the Bayesian method. The Bayesian method makes the prediction band for IG statistically meaningless. As we are interested only in the PG prediction, this is a minor problem. The IG prediction band will not be shown in the remaining results.

The prediction of PG is improved and the 95% prediction band is narrower. From a visual inspection it is seen that $S_{Y S I}$ has decreased and $S_{C G M}$ has increased, as expected. The estimates, standard deviations, and correlation table are presented in the appendix for both Data₁ and Data₂, increasing the reproducibility.

It is not meaningful to statistically compare step 1 and step 2, as information criteria are defined from the pure likelihood function; therefore visual and physiological improvements are assessed.

Model Validation

As it is desired to predict PG from IG observations, predictions without knowledge of the PG or insulin observations are computed.

Data₁

Prediction using ${\hat{θ}}_{11}$ (step 3)

Figure 3 shows the 1-step-ahead prediction using ${\hat{θ}}_{11}$ of PG using only IG observations. The prediction follows the true PG dynamic well. The 95% prediction band contains all but 1 PG observation.

Figure 3.

(Overall) Predicting with ${\hat{θ}}_{11}$ using only IG observations. Green ribbon reflects the normoglycemic range. The purple observation in the top plot is an outlier not used in the Bayesian estimation nor prediction. (Top) One-step prediction of PG and 95% prediction band. The dynamics of PG are captured in the prediction and the prediction band covers all but 1 PG observation. (Bottom) One-step prediction of IG.

Prediction using ${\hat{θ}}_{21}$ (step 6)

Figure 4 shows the predictions using ${\hat{θ}}_{21}$ and only IG observations. The dynamics of PG is captured well. The 95% prediction band contains all but 1 PG observation.

Figure 4.

(Overall) Prediction with ${\hat{θ}}_{21}$ using only IG observations. Green ribbon reflects the normoglycemic range. The purple observation in the top plot is an outlier not used in the Bayesian estimation nor prediction. (Top) One-step prediction of PG and 95% prediction band. The dynamics of PG are captured in the prediction and the prediction band covers all but 2 PG observations. (Bottom) One-step prediction of IG.

Data₂

Predicting using ${\hat{θ}}_{22}$ (step 3). Figure 5 shows Data₂, with virtually no dynamic in the PG observations and the IG observations are very fluctuating and unphysiological. Despite this difference in observations, when predicting using ${\hat{θ}}_{22},$ it is seen that the 1-step-ahead prediction of PG using only IG observations follows the PG dynamic well. The 95% prediction band contains all but 1 PG observation.

Figure 5.

(Overall) Prediction with ${\hat{θ}}_{22}$ (open-loop study) using only IG observations. Green ribbon reflects the normoglycemic range. (Top) One-step prediction of PG and 95% prediction band. Despite the limited dynamics PG is predicted well with all but 1 PG observations within the 95% prediction band. (Middle) One-step prediction of IG. (Bottom) Insulin delivery per hour.

Predicting using ${\hat{θ}}_{12}$ (step 6)

Figure 6 shows the prediction using ${\hat{θ}}_{12}$ . It is still possible to predict the dynamics of PG using only IG observations. There is a tendency that the predictions are slightly elevated compared to PG observations, but the 95% prediction band, contains all PG observation.

Figure 6.

(Overall) Prediction with ${\hat{θ}}_{12}$ (open-loop study) using only IG observations. Green ribbon reflects the normoglycemic range. (Top) One-step prediction of PG and 95% prediction band. Despite the limited dynamics PG is predicted well with all PG observations within the 95% prediction band. (Bottom) One-step prediction of IG.

Discussion

This study shows that it is possible to predict PG from IG observations using a Bayesian method. This finding indicates that the final model can be used as a control algorithm that is implementable in a pump system.

Maximum Likelihood (Step 1)

Using maximum likelihood estimation, the IG influence on the PG prediction is evident. The many closely spaced observations have a tremendous influence on the likelihood.³² The optimization minimizes the 1-step prediction, and as there are more IG observations, the optimizer will gain more from explaining the IG dynamics. To avoid this, a Bayesian method is suggested as a solution.

A Bayesian Method (Steps 2, 5 and 6)

It is assumed that the global parameter set for a patient will vary between days. This variation is assumed to have a limited span. This assumption enables us to cross-validate within the same patient. We expect that the parameter estimates from a given night will be similar to parameter estimates from another night. This allows us to use prior knowledge from 1 study night to estimate a new parameter set.

The ideal solution for a Bayesian method implementation is estimation of the optimal prior knowledge. However, an estimation of the variance of the variance of the measurement device requires a large amount of data. Instead, the literature was used to approximate measurement-device prior knowledge. The 95% prediction band is wide on the IG prediction since a Bayesian method is used. But as the focus of the study is prediction of PG this is accepted. The 95% prediction bands for PG are narrow.

As discussed by Donnet and Samson, Bayesian methods have been proposed and implemented for population pharmacokinetics/pharmacodynamics SDE-GB.⁴³ Using CTSM-R it is possible to implement Bayesian estimation on a single subject, obtaining valid and interpretable models. With this approach insulin observations are needed to estimate the parameters of the model. A similar approach could have been tested using clinical knowledge to fix the parameter. Future studies are needed to clarify this and should include several patients and preferably population modeling.

Physiological Modeling

The SDE-GB models applied allow for a physiological, opposed to a purely statistical, interpretation of the models, since the model structure is based on physiology. From a visual inspection it is evident that the statistical best model is not the physiologically and technologically most correct model.

The use of CTSM-R allows optimal conditions for the modeling procedure: that is, continuous-time dynamics with discrete time observations and a possibility of adopting a Bayesian modeling framework. This gives fewer parameters to estimate, which in turn gives a more robust estimation. Furthermore, previous knowledge is often given as differential equations, and therefore easily incorporated into existing models.⁴³ This maintains the physiological meaning of the model, and gives it interpretability for both the modelers and medical experts.

The examination of the covariance matrixes reveals that there is no severe parameter correlation. The covariance matrix has great potential in constructing a virtual patient database. Randomly selected parameter sets can be used to construct an infinite number of simulations for use as initial investigations, instead of costly human experiments.

Perspectives

A natural progression is to extend the Bayesian method to data where other factors, such as meals, have an influence. Building on Reifman et al.,¹⁶ a combination of data-driven AR models and Bayesian probability models should be investigated. The Bayesian probability ensures physiological interpretability and the AR processes capture sensor dynamics.

New optical-sensor techniques are currently under development as an alternative to the electrochemical transducer. Vaddiraju et al. provide an overview of advantages, disadvantages, and perspectives.⁷

To obtain more knowledge about PG-IG dynamics, experiments with multiple glucose sensors, both electrochemical and optical, should be considered. Advantages of a multiple glucose sensor approach are not well investigated, but reports describing the advantages exist.^44,45

Conclusion

The aim of this article was to extend an existing model describing PG-IG dynamical behavior to SDE-GB form With the SDE-GB model, and clinical data from a closed-loop study, parameters were estimated. The initial investigation used a maximum likelihood approach. Afterward a Bayesian method was investigated.

Using the maximum likelihood approach gave statistically reliable result, but there were deficiencies, as the ability to predict PG from IG observations was poor. A Bayesian method was successfully implemented, and it has been shown that parameters can be estimated using CTSM-R. Using a Bayesian estimate we were able to predict PG using IG observations. Furthermore we were able to cross-validate parameter estimates using a different data set from the same patient. Using prior knowledge from a previous study night, we were able to estimate usable parameters, indicating a close relation in intrapatient parameters.

These findings indicate that the model can be incorporated in an artificial pancreas that can be used as a closed-loop controller.

Footnotes

Appendix.

Estimated Parameter Values.

Data₁						Correlation table
Parameters	Unit	Prior $p r i o r (θ_{M D})$	Estimate	Lower limit 2.5 %	Upper limit 97.5 %	EGP	GEZI	p2	$e_{C G M}$	Si	$σ_{I_{S C}}$	$σ_{G_{P}}$	$σ_{G_{S C}}$	$e_{I N S}$	$e_{Y S I}$	$τ_{1}$	$τ_{3}$
EGP	mg/(dL×min)		0.0504*	0.0037	0.075	1
GEZI	1/min		0.0000001	0	0.00036	−.0002	1
p2	1/min		0.0098*	0.00069	0.0137	−.5797	.017	1
$S_{C G M}$	—	1.04971	0.4481*	0.408	0.492	.0033	−.005	.000	1
Si	L/(mU×min)		0.3772*	0.245	0.582	.9509	−.008	−.709	.003	1
$σ_{I_{S C}}$	—		0.0012*	0.000070	0.0020	.0063	−.004	.028	−.006	−.004	1
$σ_{G_{P}}$	—		0.0636*	0.00449	0.089	.1266	−.020	−.064	−.003	.132	−.008	1
$σ_{G_{S C}}$	—		0.1040*	0.00731	0.1479	−.0231	.022	−.010	.019	−.021	.046	.006	1
$S_{I N S}$	—	20.0855	0.0000*	0.000	0.000	.0076	−.009	−.040	−.010	.020	−.054	−.019	.057	1
$S_{Y S I}$	—	1.12277	0.0048*	0.000387	0.0059	−.0021	−.033	.016	−.011	.001	.018	−.005	.018	−.005	1
$τ_{1}$	min		65.5687*	54.09	75.90	−.0182	.026	.025	−.006	−.009	−.196	−.005	.039	.262	.032	1
$τ_{3}$	min		42.6466*	9.46	75.53	.0500	−.030	.016	−.042	.047	−.008	−.035	−.337	−.007	−.092	.021	1
Data₂						Correlation table	Data₂						Correlation table	Data₂
Parameters	Unit	Prior $p r i o r (θ_{M D})$	Estimate	Lower limit 2.5 %	Upper limit 97.5 %	EGP	Parameters	Unit	Prior $p r i o r (θ_{M D})$	Estimate	Lower limit 2.5 %	Upper limit 97.5 %	EGP	Parameters	Unit	Prior $p r i o r (θ_{M D})$	Estimate
EGP	mg/(dL×min)		0.0392*	0.0199	0.0769	1	EGP	mg/(dL×min)		0.0392*	0.0199	0.0769	1	EGP	mg/(dL×min)		0.0392*
GEZI	1/min		0.0091*	0.0034	0.0242	.933	GEZI	1/min		0.0091*	0.0034	0.0242	.933	GEZI	1/min		0.0091*
p2	1/min		0.0573	0.0019	1.6597	−.735	p2	1/min		0.0573	0.0019	1.6597	−.735	p2	1/min		0.0573
$S_{C G M}$	—	1.04971	0.0503*	0.0015	1.6437	.005	$S_{C G M}$	—	1.04971	0.0503*	0.0015	1.6437	.005	$S_{C G M}$	—	1.04971	0.0503*
Si	L/(mU×min)		0.4208	0.3854	0.4594	−.549	Si	L/(mU×min)		0.4208	0.3854	0.4594	−.549	Si	L/(mU×min)		0.4208
$σ_{I_{S C}}$	—		0.0022*	0.0014	0.0033	−.003	$σ_{I_{S C}}$	—		0.0022*	0.0014	0.0033	−.003	$σ_{G_{P}}$	—		0.0022*
$σ_{G_{P}}$	—		0.0247*	0.0160	0.0380	.157	$σ_{I_{S C}}$	—		0.0247*	0.0160	0.0380	.157	$σ_{G_{S C}}$	—		0.0247*
$σ_{G_{S C}}$	—		0.0975*	0.0589	0.1613	−.005	$σ_{G_{S C}}$	—		0.0975*	0.0589	0.1613	−.005	$σ_{G_{P}}$	—		0.0975*
$S_{I N S}$	—	20.0855	0.0000*	0.000	0.000	.022	$S_{I N S}$	—	20.0855	0.0000*	0.000	0.000	.022	$S_{I N S}$	—	20.0855	0.0000*
$S_{Y S I}$	—	1.12277	0.0047*	0.0038	0.0058	.011	$S_{Y S I}$	—	1.12277	0.0047*	0.0038	0.0058	.011	$S_{Y S I}$	—	1.12277	0.0047*
$τ_{1}$	min		54.0378*	28.52	79.48	.006	$τ_{1}$	min		54.0378*	28.52	79.48	.006	$τ_{1}$	min		54.0378*
$τ_{3}$	min		74.2980*	11.28	136.72	−.006	$τ_{3}$	min		74.2980*	11.28	136.72	−.006	$τ_{3}$	min		74.2980*

The major differences in prior knowledge come from a transformation to log scale in CTSM-R where the value of the prior is based on the initial value of the estimation.

Significant at α = .05.

Abbreviations

CGM, continuous glucose monitor; CSII, continuous subcutaneous insulin infusion; IG, interstitial glucose; MVP, Medtronic virtual patient; PG, plasma glucose; SDE-GB, stochastic-differential-equation-based gray-box; SDEs, stochastic differential equations; T1DM, type 1 diabetes mellitus.

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: This work was funded as a part of the DIACON project by the Danish Council of Strategic Research (NABIIT project 2106-07-0034).

References

Pickup

Freeman

Sutton

. Glycaemic control in type 1 diabetes during real time continuous glucose monitoring compared with self monitoring of blood glucose: meta-analysis of randomised controlled trials using individual patient data. BMJ. 2011;343:d3805-d3805.

Bequette

. Continuous glucose monitoring: real-time algorithms for calibration, filtering, and alarms. J Diabetes Sci Technol. 2010;4(2):404-418.

Keenan

Mastrototaro

Voskanyan

Steil

. Delays in minimally invasive continuous glucose monitoring devices: a review of current technology. J Diabetes Sci Technol. 2009;3(5):1207-1214.

Kulcu

Tamada

Reach

Potts

Lesho

. Physiological differences between interstitial glucose and blood glucose measured in human subjects. Diabetes Care. 2003;26(8):2405-2409.

Steil

Rebrin

Mastrototaro

Bernaba

Saad

. Determination of plasma glucose during rapid glucose excursions with a subcutaneous glucose sensor. Diabetes Technol Ther. 2003;5(1):27-31.

Rebrin

Sheppard

Jr Steil

. Use of subcutaneous interstitial fluid glucose to estimate blood glucose: revisiting delay and sensor offset. J Diabetes Sci Technol. 2010;4(5):1087-1098.

Vaddiraju

Burgess

Tomazos

Jain

Papadimitrakopoulos

. Technologies for continuous glucose monitoring: current problems and future promises. J Diabetes Sci Technol. 2010;4(6):1540-1562.

Kumareswaran

Elleri

Allen

. Accuracy of continuous glucose monitoring during exercise in type 1 diabetes pregnancy. Diabetes Technol Ther. 2013;15(3):223-229.

Gerritsen

Jansen

Lutterman

. Performance of subcutaneously implanted glucose sensors for continuous monitoring. Neth J Med. 1999;54(4):167-179.

10.

Wilson

Beck

Tamborlane

. The accuracy of the FreeStyle Navigator continuous glucose monitoring system in children with type 1 diabetes. Diabetes Care. 2007;30(1):59-64.

11.

Diabetes Research in Children Network (DirecNet) Study Group. Lack of accuracy of continuous glucose sensors in healthy, nondiabetic children: results of the Diabetes Research in Children Network (DirecNet) accuracy study. J Pediatr. 2004;144(6):770-775.

12.

Phillip

Battelino

Atlas

. Nocturnal glucose control with an artificial pancreas at a diabetes camp. N Engl J Med. 2013;368(9):824-833.

13.

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.

14.

Kovatchev

Clarke

. Peculiarities of the continuous glucose monitoring data stream and their impact on developing closed-loop control technology. J Diabetes Sci Technol. 2008;2(1):158-163.

15.

Hoss

Jeddi

Schulz

Budiman

Bhogal

McGarraugh

. Continuous glucose monitoring in subcutaneous tissue using factory-calibrated sensors: a pilot study. Diabetes Technol Ther. 2010;12(8):591-597.

16.

Reifman

Rajaraman

Gribok

Ward

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

17.

Knobbe

Buckingham

. The extended Kalman filter for continuous glucose monitoring. Diabetes Technol Ther. 2005;7(1):15-27.

18.

Facchinetti

Sparacino

Cobelli

. Enhanced accuracy of continuous glucose monitoring by online extended Kalman filtering. Diabetes Technol Ther. 2010;12(5):353-363.

19.

Leal

Garcia-Gabin

Bondia

. Real-time glucose estimation algorithm for continuous glucose monitoring using autoregressive models. J Diabetes Sci Technol. 2010;4(2):391-403.

20.

Chase

Hann

Jackson

. Integral-based filtering of continuous glucose sensor measurements for glycaemic control in critical care. Comput Methods Programs Biomed. 2006;82(3):238-247.

21.

Breton

Kovatchev

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

22.

Guerra

Facchinetti

Sparacino

De Nicolao

Cobelli

. Enhancing the accuracy of subcutaneous glucose sensors: a real-time deconvolution-based approach. IEEE Trans Biomed Eng. 2012;59(6):1658-1669.

23.

Pickup

Williams

. Textbook of Diabetes. 2 vols. New York, NY: John Wiley; 2003.

24.

Madsen

. Time Series Analysis. Boca Raton, FL: CRC Press; 2008.

25.

Møller

Overgaard

Madsen

Hansen

Pedersen

Ingwersen

. Predictive performance for population models using stochastic differential equations applied on data from an oral glucose tolerance test. J Pharmacokinet Pharmacodyn. 2009;37(1):85-98.

26.

Duun-Henriksen

Schmidt

Røge

. Model identification using stochastic differential equation grey-box models in diabetes. J Diabetes Sci Technol. 2013;7(2):431-440.

27.

Kristensen

Madsen

Jørgensen

. A method for systematic improvement of stochastic grey-box models. Comput Chem Eng. 2004;28(8):1431-1449.

28.

Tornøe

Overgaard

Agersø

Nielsen

Madsen

Jonsson

. Stochastic differential equations in NONMEM®: implementation, application, and comparison with ordinary differential equations. Pharm Res. 2005;22(8):1247-1258.

29.

Overgaard

Jonsson

Tornøe

Madsen

. Non-linear mixed-effects models with stochastic differential equations: implementation of an estimation algorithm. J Pharmacokinet Pharmacodyn. 2005;32(1):85-107.

30.

Kristensen

Madsen

Ingwersen

. Using stochastic differential equations for PK/PD model development. J Pharmacokinet Pharmacodyn. 2005;32(1):109-141.

31.

Björk

. Arbitrage Theory in Continuous Time. Oxford, UK: Oxford University Press; 2004.

32.

Madsen

. Introduction to General and Generalized Linear Models. Boca Raton, FL: CRC Press; 2011.

33.

Kanderian

Weinzimer

Voskanyan

Steil

. Identification of intraday metabolic profiles during closed-loop glucose control in individuals with type 1 diabetes. J Diabetes Sci Technol. 2009;3(5):1047-1057.

34.

Schmidt

Boiroux

Duun-Henriksen

. Model-based closed-loop glucose control in type 1 diabetes—the DiaCon experience. J Diabetes Sci Technol. 2013;7(5):1255-1264.

35.

Bergman

Phillips

Cobelli

. Physiologic evaluation of factors controlling glucose tolerance in man: measurement of insulin sensitivity and beta-cell glucose sensitivity from the response to intravenous glucose. J Clin Invest. 1981;68(6):1456-1467.

36.

Lynch

Bequette

. Model predictive control of blood glucose in type I diabetics using subcutaneous glucose measurements. Am Control Conf 2002 Proc. 2002;5:4039-4043.

37.

Fisher

. A semiclosed-loop algorithm for the control of blood glucose levels in diabetics. IEEE Trans Biomed Eng. 1991;38(1):57-61.

38.

Wilinska

Chassin

Schaller

Schaupp

Pieber

Hovorka

. Insulin kinetics in type-1 diabetes: continuous and bolus delivery of rapid acting insulin. IEEE Trans Biomed Eng. 2005;52(1):3-12.

39.

Fisher

. On the mathematical foundations of theoretical statistics. Philos Trans R Soc Lond Ser Contain Pap Math Phys Character. 1922;222(594-604):309-368.

40.

Khan

Vasquez

Gray

Wians

Jr Kroll

. The variability of results between point-of-care testing glucose meters and the central laboratory analyzer. Arch Pathol Lab Med. 2006;130(10):1527-1532.

41.

Dexcom SEVEN PLUS user guide. http://www.dexcom.com/sites/dexcom.com/files/seven-Plus/docs/seven_plus_users_guide.pdf.

42.

CTSMR. www.ctsm.info.

43.

Donnet

Samson

. A review on estimation of stochastic differential equations for pharmacokinetic/pharmacodynamic models. Adv Drug Deliv Rev. Available at: http://www.sciencedirect.com/science/article/pii/S0169409X13000501. Accessed June 17, 2013.

44.

Harman-Boehm

Gal

Raykhman

Naidis

Mayzel

. Noninvasive glucose monitoring: increasing accuracy by combination of multi-technology and multi-sensors. J Diabetes Sci Technol. 2010;4(3):583-595.

45.

Castle

Ward

. Amperometric glucose sensors: sources of error and potential benefit of redundancy. J Diabetes Sci Technol. 2010;4(1):221-225.

Predicting Plasma Glucose From Interstitial Glucose Observations Using Bayesian Methods

Abstract

Background:

Method:

Results:

Conclusions:

Keywords

Methods

The Initial Model

The Maximum Likelihood Principle

Modeling Process

Step 1 (Diffusion Term Identification and Maximum Likelihood)

Steps 2-6 (A Bayesian Method)

Step 2

Step 3

Step 4

Step 5

Step 6

Results

Maximum Likelihood (Step 1)

A Bayesian Method (Step 2)

Model Validation

Data1

Prediction using θ ^ 11 (step 3)

Prediction using θ ^ 21 (step 6)

Data2

Predicting using θ ^ 12 (step 6)

Discussion

Maximum Likelihood (Step 1)

A Bayesian Method (Steps 2, 5 and 6)

Physiological Modeling

Perspectives

Conclusion

Footnotes

Appendix.

Abbreviations

Declaration of Conflicting Interests

Funding

References

Data₁

Prediction using ${\hat{θ}}_{11}$ (step 3)

Prediction using ${\hat{θ}}_{21}$ (step 6)

Data₂

Predicting using ${\hat{θ}}_{12}$ (step 6)