Plasma Insulin Estimation in People with Type 1 Diabetes Mellitus

Article pubs.acs.org/IECR

Plasma Insulin Estimation in People with Type 1 Diabetes Mellitus Iman Hajizadeh,† Mudassir Rashid,† Kamuran Turksoy,‡ Sediqeh Samadi,† Jianyuan Feng,† Nicole Frantz,‡ Mert Sevil,‡ Eda Cengiz,§ and Ali Cinar*,† †

Department of Chemical and Biological Engineering and ‡Department of Biomedical Engineering, Illinois Institute of Technology, Chicago, Illinois 60616, United States § Department of Pediatrics, Yale University School of Medicine, New Haven, Connecticut 06437-2411, United States S Supporting Information *

ABSTRACT: In this work the real-time estimation of plasma insulin concentration (PIC) to quantify the insulin in the bloodstream in patients with type 1 diabetes mellitus (T1DM) is presented. To this end, Hovorka’s model, a glucose−insulin dynamics model, is incorporated with various estimation techniques, including continuous-discrete extended Kalman filtering, unscented Kalman filtering, and moving horizon estimation, to provide an estimate of PIC. Furthermore, due to the considerable variability in the temporal dynamics of patients, some uncertain model parameters that have significant effects on PIC estimates are considered as additional states in Hovorka’s model to be simultaneously estimated. Latent variable regression models are developed to individualize the PIC estimators by appropriately initializing the time-varying model parameters for improved convergence. The performance of the proposed methods is evaluated using clinical data sets from subjects with T1DM, and the results demonstrate the accurate estimation of PIC.

■

INTRODUCTION In healthy people, the pancreas produces insulin that regulates blood glucose levels. In people with type 1 diabetes mellitus (T1DM), insulin production ceases because an autoimmune attack destroys the insulin-producing cells (β-cells) in the pancreas. Although people suffering from T1DM routinely take insulin injections to regulate their blood glucose concentration (BGC), they still face numerous complications due to high BGC (i.e., damage to blood vessel, kidneys, nerve, and eye) or low BGC (i.e., dizziness, unconsciousness, seizures, and coma). Closed-loop control of BGC, known as the artificial pancreas (AP) system, is an alternative treatment for people with T1DM to automatically regulate their blood glucose levels.1−8 A basic AP system uses continuous glucose monitoring (CGM) measurements to calculate the optimum amount of insulin to be infused with an insulin pump in order to maintain glucose levels within a safe target range. One of the common methods that an AP system uses to control the blood glucose levels is basal−bolus insulin therapy. Basal insulin is a constant low dose of insulin to keep blood glucose levels at consistent levels during periods of fasting. Bolus insulin is a single large dose of insulin that is usually administered at meal times to control blood glucose levels after a meal. The injected insulin (basal or bolus) gradually accumulates in the bloodstream and is eventually utilized by the body. One of the factors that prolongs the utilization of administered insulin is the significant time delays involved in the diffusion and absorption of the subcutaneously injected insulin analogues. Furthermore, increased insulin resistance (typically caused by © 2017 American Chemical Society

low physical activity levels) is often accompanied by elevated concentrations of plasma insulin. Hence, an AP control system that does not explicitly consider these important factors may administer high amounts of insulin before the previously administered insulin present in the bloodstream is completely utilized and its effects are observed in the CGM sensor measurements. Consequently, fear of low BGC due to the overdosing of insulin is a major concern for people with T1DM using an AP system, and a few closed-loop studies involving various control algorithms have resulted in mild or severe hypoglycemic episodes.9 Hence, it is crucial for an AP system to be cognizant of the amount of insulin already present in the bloodstream and to adjust the future insulin doses accordingly in order to avoid episodes of hypoglycemia. Despite their necessity, sensors that can accurately measure the amount of insulin in the blood, or plasma insulin concentration (PIC), in real time are not readily available. Some common methods for incorporating bloodstream insulin information in various AP control strategies involve consulting insulin-on-board (IOB) curves, estimating plasma insulin concentrations, and calculating the amount of subcutaneously administered insulin through insulin absorption models.1,2,10−16 One widely employed approach is to consider the amount of active insulin in the body from previous doses, Received: Revised: Accepted: Published: 9846

April 17, 2017 July 7, 2017 August 16, 2017 August 16, 2017 DOI: 10.1021/acs.iecr.7b01618 Ind. Eng. Chem. Res. 2017, 56, 9846−9857

Article

Industrial & Engineering Chemistry Research

are designed based on Hovorka’s glucose−insulin model.25 The main challenges in the design of PIC estimators involving Hovorka’s glucose−insulin model and using CGM and infused insulin data that need to be addressed are (i) the observability of model states, (ii) the discrete sampling of the measurement output (CGM), (iii) the limited knowledge about the time and exact amount of meals and other disturbances, (iv) the timevarying nature of human physiology, (v) the intersubject variability, and (vi) the constraints on the state variables and the rate of change in the parameters of the model. In this study, three different methods are evaluated to simultaneously estimate the states and parameters: continuousdiscrete extended Kalman filter (CDEKF), unscented Kalman filter (UKF), and moving horizon estimator (MHE). In addition, the effect of meal intake UG is considered a time-varying model parameter to be estimated by these techniques as well. Therefore, the proposed algorithms are able to estimate PIC without requiring onerous and obscure information, such as the exact time or amount of meal intake, by relying on the real-time CGM measurements and infused insulin data. Because of the significant temporal variability in the glucose−insulin dynamics, two unknown model parameters that have a significant effect on PIC, insulin elimination from plasma (ke) and time-to-maximum of the absorption of subcutaneously injected insulin (tmax,I), are also defined to be time-varying parameters to be estimated online so that the intrapatient variability is well characterized.17 Hence, the parameters tmax,I, ke, and UG (listed along other parameters in Table 1) are considered as extended states in Hovorka’s glucose−insulin model for the CDEKF and UKF algorithms and manipulated variables in the MHE formulation. Providing appropriate initial conditions for the parameters tmax,I and ke due to the intersubject variability of people with T1DM is necessary to obtain rapid convergence of the estimated parameters online. To determine appropriate initial values for these unknown model parameters, partial least-squares (PLS) models based on each patient’s demographic information such as body weight, height, BMI, and total daily insulin dose are developed. In order to incorporate constraints on the estimated states (such as bounds to avoid physically unrealizable or unrealistic values) and to restrict the maximum allowable change in the time-varying parameters at each sampling time, modified versions of the estimation algorithms are utilized. The paper is organized as follows. The glucose−insulin dynamics model (Hovorka’s model) with all variables, parameters, units, and nominal values is presented in the next section, followed by a brief review of observability. Then three different techniques, namely, the CDEKF, UKF, and MHE algorithms for simultaneous state and parameter estimation, and partial least-squares models for the individualization of PIC estimators are presented along with detailed information about the subjects and study experiments. Finally, results, discussions, and conclusions are provided.

referred to as insulin-on-board, via predetermined IOB curves that represent the time duration of insulin action in the body. Consulting these curves assists patients in discontinuing the infusion of additional insulin if there is sufficient insulin already present in the bloodstream, thus avoiding insulin stacking that could result in the overcorrection of the BGCs and may lead to hypoglycemia. These curves have also been used in different AP control strategies to limit of insulin infusion. Nevertheless, this information from IOB curves is an approximate prediction of active insulin in the body and is not a direct estimate of plasma insulin concentration. Although IOB curves provide valuable information about the active amount of insulin in the body, the significant time-varying delays induced by the subcutaneous administration of insulin and changes in a person’s metabolic state throughout the day affect insulin utilization and the IOB. Therefore, the accurate estimation of PIC in real time can provide important information to AP systems for calculating the optimum amount of insulin to be infused. Furthermore, estimates of PIC can be obtained from adaptive observers designed for simultaneous state and parameter estimation based on reliable glucose−insulin models and using real measured outputs (CGM data) from patients. Several methods have been proposed for real-time estimation of PICs.8,17−24 An extended Kalman filter was designed to provide real-time estimates of PIC from CGM data based on Hovorka’s glucose−insulin model with various time-varying model parameters considered as extended states in Hovorka’s original model.17 The proposed method was tested using an insilico study of 100 patients with T1DM and clinical data from 12 insulin pump patients demonstrating satisfactory estimation results. An unscented Kalman filter was designed to provide an estimate of the current PIC based on the measurement of the plasma glucose using a discretized version of Bergman’s minimal model with plasma glucose concentration as an output.18 Data from an intravenous glucose tolerance test (IVGTT) was used to evaluate the estimation results. Moreover, the unscented Kalman filter approach was also applied on the extended Bergman’s model to simultaneously estimate parameters and states.19 An estimator incorporating error feedback was also proposed based on the measured and predicted plasma glucose using Bergman’s third-order nonlinear model designed to tolerate measurement noise as well as discretization errors by means of the H∞ criterion.20 The proposed H∞ estimator was tested using synthetic patient data produced using IVGTT. A particle filter approach using Bergman’s minimal model is proposed to estimate PIC, and the data set from an IVGTT was used to evaluate the performance.21 Additionally, a nonlinear observer was designed on the basis of an aggregate model where Bergman’s nonlinear model was expanded in terms of local linear models weighted by normalized radial basis functions.22 The performance of the nonlinear observer is assessed using the data set from IVGTT as well. This work considers the problem of designing reliable PIC estimators that take into account the significant temporal variability in the glucose−insulin dynamics. This variability manifests itself in both inter- and intrasubject fluctuations and is an important consideration with regard to the performance of AP systems. Individuals with T1DM exhibit a diverse range of responses to carbohydrate intake and insulin infusion, with significant differences among patients and also considerable variations within individuals. Thus, the control law in an AP system requires a reliable, time-varying, and individualized PIC estimator to perform efficiently. In this study, the PIC estimators

■

METHODS Glucose−Insulin Dynamic Model. Several models of insulin action and glucose kinetics using clinical data have been proposed to describe glucose production and utilization and insulin and meal absorption. The Dalla Man-Cobelli model,26 Hovorka’s model,25 and the Bergman minimal model27 are some of the widely used physiological models of insulin action and the glucose kinetics system for developing various control or estimation techniques. In this study, Hovorka’s model, detailed in eq 1 for completeness, is used for designing the estimators. 9847

DOI: 10.1021/acs.iecr.7b01618 Ind. Eng. Chem. Res. 2017, 56, 9846−9857

Article


equation. This model contains nine states as described by the system given in eq 1.

Table 1. Variables, Parameters, and Nominal Values for the Hovorka’s Model25 symbol

definition

S1,S2

Gsub UG u k12 ka1 ka2 ka3 ke VI VG SfIT

subcutaneously administered insulin absorption plasma insulin concentration effects of insulin on glucose distribution/ transport effects of insulin on glucose disposal effects of insulin on endogenous glucose production masses of glucose in the accessible compartment masses of glucose in the nonaccessible compartment subcutaneous glucose gut absorption rate administration (bolus and basal) of insulin transfer rate deactivation rate deactivation rate deactivation rate insulin elimination from plasma insulin distribution volume glucose distribution volume insulin sensitivity of distribution/transport

SfID

insulin sensitivity of disposal

SfIE EGP0

insulin sensitivity of EGP EGP extrapolated to zero insulin concentration noninsulin-dependent glucose flux

I x1 x2 x3 Q1 Q2

F01 tmax,I BW kb1 kb2 kb3 Fc01 FR

dS1(t ) S (t ) = u(t ) − 1 dt tmax, I

value (units) − (mU) − (mU L−1) − (min−1)

dS2(t ) S (t ) S (t ) = 1 − 2 dt tmax, I tmax, I

− (min−1) − (min−1)

S (t ) d I (t ) = 2 − keI(t ) dt tmax, IVI

− (mmol)

dx1(t ) = −ka ,1x1(t ) + kb ,1I(t ) dt

− (mmol) − (mmol L−1) − (mmol min−1) − (mU min−1) 0.066 (min−1) 0.006 (min−1) 0.06 (min−1) 0.03 (min−1) 0.138 (min−1) 0.12 × BW (L) 0.16 × BW (L) 51.2 × 10−4(L min−1 mU−1) 8.2 × 10−4 (L min−1 mU−1) 520 × 10−4 (L mU−1) 0.0161 (mmol kg−1 min−1) 0.0097 (mmol kg−1 min−1) 55 (min)

dx 2(t ) = −ka ,2x 2(t ) + kb ,2I(t ) dt dx3(t ) = −ka ,3x3(t ) + kb ,3I(t ) dt dQ 1(t ) dt

c (t ) + FR(t ) − x1(t )Q 1(t ) = UG(t ) − F01

+ k12Q 2(t ) + EGP0(1 − x3(t )) dQ 2(t ) dt

= x1(t )Q 1(t ) − (k12 + x 2(t ))Q 2(t )

⎞ dGsub(t ) 1 ⎛ Q (t ) = ⎜ 1 − Gsub(t )⎟ dt τ ⎝ VG ⎠

(1)

Plasma Insulin Estimator. The dynamic model (eq 1) can be expressed in the following form

time-to-maximum of absorption of injected insulin weight − (kg) =SfIT × ka1 f =SID × ka2 f =SIE × ka3 F01 if G ≥ 4.5 mmol/L; otherwise F01G/4.5 0.003 (G−9)VG if G ≥ 9 mmol/L; otherwise 0

d X (t ) = f (X(t ), u(t )) dt y(t ) = h(X(t ))

(2)

where X(t) = [S1(t), S2(t), I(t), x1(t), x2(t), x3(t), Q1(t), Q2 (t), Gsub(t)] denotes the vector of state variables and h(X(t)) denotes the measurement function with y(t) as the output of the model, where the output measurement variable is the subcutaneous glucose Gsub(t). When considering the subsystem that describes the subcutaneous insulin infusion and PIC, including the dynamics of the state variables S1(t), S2(t), and I(t), the parameters tmax,I and ke have a direct effect on the PIC. In this study, ke and tmax,I are uncertain parameters and considered as extended states in model 2. Furthermore, as information about the time and quantity of meals is difficult to ascertain and thus considered as unknown disturbances, the gut absorption rate UG(t) is also included as an extended state in model 2. Incorporating the uncertainty in the dynamics of the system (referred to as process noise) and the noise in the discrete-time measurements28 we have

This model has differential equations describing the glucose− insulin dynamics, the subsystem describing the blood glucose dynamics, the subsystem describing the subcutaneous insulin infusion, and the subsystem modeling the glucose transport from plasma to interstitial tissues. All variables, parameters, units, and nominal values used in this model (eq 1) are introduced in Table 1.25 The blood glucose concentration dynamics are described using a two-compartment model. The two state variables Q1(t) and Q2(t) denote the masses of glucose in the accessible and nonaccessible compartments, respectively. Another two-compartment model, with state variables S1(t) and S2(t), defines the absorption rate of subcutaneously administered insulin. The PIC, I(t), is represented by a first-order differential equation. The insulin action is calculated by using three variables, the influence on transport and distribution x1(t), the utilization and phosphorylation of glucose in adipose tissue x2(t), and the endogenous glucose production in the liver x3(t). The relationship between BGC and subcutaneous glucose reported by CGM measurements is considered a first-order dynamic

dX ′(t ) = f ′(X ′(t ), u(t )) + G(t )ω(t ), ω(t ) ≈ N (0, Q ) dt yk = h(Xk′) + νk , νk ≈ N (0, R ) (3)

where ω(t) and νk represent the process and observation noise vectors, Q(t) and R(t) denote the covariance and variance of these noises, respectively, and X′(t) = [S1(t), S2(t), I(t), x1(t), x2(t), x3(t), Q1(t), Q2 (t), Gsub(t), tmax,I (t), ke(t), UG(t)] denotes 9848


Article

Industrial & Engineering Chemistry Research the augmented state vector including the uncertain model parameters. To construct the estimator (observer), the system described by eq 3 should be observable. Methods to check the observability for nonlinear systems include using empirical observability Gramian or calculating the observability matrix using the Lie derivatives of the output function. The observability of the original Hovorka’s model and with extended states for real-time estimation of tmax,I and ke have already been proven.17 In this study, the observability of the extended Hovorka’s model with the augmented state vector including the gut absorption rate UG(t) is investigated. The following observability matrix is used to prove the observability of the model eq 3 ⎡ ∂ 0 ⎤ Lf ′(h(X ′(t ))) ⎥ ⎢ ∂ X ′ ⎢ ⎥ ⎢ ∂ 1 ⎥ WO(t ) = ⎢ ∂X ′ Lf ′(h(X ′(t ))) ⎥ ⎢ ⎥ ... ⎢ ⎥ ⎢ ∂ n−1 ⎥ Lf ′ (h(X ′(t )))⎥ ⎢⎣ ⎦ ∂X ′

Table 2. Continuous-Discrete Extended Kalman Filter initialization time update

measurement update

Jacobians

P(0) = P0, X̂ ′(0) = X̅ 0, where P is covariance matrix estimate: X̂ k′− from model Ẋ̂ ′(t) = f(X̂ ′(t),u(t)) error covariance: Ṗ(t) = A(X̂ ′(t))P(t) + P(t)A(X̂ ′(t))T + GQ(t)G Kalman gain: Kk = P−k HT(X̂ k′−)[H(X̂ k′−)P−k HT(X̂ k′−) + R]−1 error covariance: Pk = [I − KkH(X̂ k′−)]P−k estimate: X̂ ′k = X̂ ′k− + Kk(yk − h(X̂ ′k−)) A(X̂ ′(t))=

∂f (X̂ ′ (t )) , ∂X ′

H(X̂ k′)=

∂h(Xk̂ ′) ∂X ′

reasonable bounds that ensure physically realizable values are obtained. For example, an estimator designed by Hovorka’s model mathematically may report a negative value for PIC, but this is not physically possible. Hence, the use of constraints is important to ensure the estimates correspond with the physical definitions. One of the widely used methods for implementing constraints in the EKF technique is known as clipping where the corrected estimator states are set equal to predefined upper or lower bounds if they are outside those specified limits.30−35 In this study, the clipping technique is used to handle the constrains in CDEKF. In addition to the bounds on the estimated states, maximum rates of change constraints are defined for the parameters tmax,I and ke. The reason for the rate of change constraint is to avoid sudden changes in these parameters due to high noise levels or unknown disturbances at certain sampling times that may result in inappropriate corrections. All imposed constraints, including the bounds of the state estimates and the rate of change constraints, are presented in eq 5

(4)

where n = 12 is the order of system 3, Lkf′(h(X′(t))) ≔ ∂Lkf ′− 1(h(X ′ (t )))

f ′, and L0f′(h(X′(t))) ≔ h(X′(t)). The observability matrix 4 was calculated, and it resulted in being of rank n, so the system described by eq 3 is observable as well. The rewritten model 3 is used by the CDEKF and UKF methods. For developing the MHE formulation, the original Hovorka’s model described by eq 1 is used with the uncertain parameters considered to be additional decision variables in the optimization. Continuous-Discrete Extended Kalman Filter. The extended Kalman filter was originally developed for discrete-time systems, yet most physical systems are represented with continuous-time models, while the outputs are sampled in a discrete manner. Hence, the formulation of continuous-discrete extended Kalman filtering is appealing for such systems where the state evolution model is continuous with discrete output sampling. Since Hovorka’s glucose−insulin model is a nonlinear continuous model and the information for CGM measurements are available every 5 min, a continuous-discrete extended Kalman filter is appropriate for this application. This type of Kalman filter includes two steps, prediction and correction based on discretetime measurements, and it is able to provide estimation for the unmeasured states continuously. Using a continuous-time insulin−glucose model provides the advantage of having information available at higher sampling rates for use in various AP modules such as performance monitoring, fault detection and diagnosis, and control system performance assessment modules that enable the AP to perform more efficiently. The other advantage of this method is the elimination of the model discretization step, thus avoiding the error resulting from the discrete approximation of the continuous-time model that might negatively influence the accuracy of the estimator. In comparison to the other estimation methods such as UKF and MHE, this method is faster in real time since the amount of necessary computations is less than the other approaches. The CDEKF is designed by using the model described by eq 3. A brief introduction of this technique is presented in Table 2.29 As the states represent a physical process through firstprinciples models, the states should be maintained within ∂X ′

Lb ≤ X̂k′ ≤ Ub k k−1 |tmax, I − tmax, I | ≤ δ1

|kek − kek − 1| ≤ δ2

(5)

Unscented Kalman Filtering. In the literature, some limitations have been mentioned for the EKF.30 For example, since the nonlinear system is linearized in the EKF algorithm for the state prediction and covariance matrix calculation, the derived linear system may not accurately describe the nonlinear system dynamics at some sampling times, which may result in instability or poor performance. The unscented Kalman filter can handle these limitations by using the unscented transformation (UT). The UT is a method for calculating the statistics of a random variable that undergoes a nonlinear transformation. The UT characterizes the mean and covariance estimates with a minimal set of sample points called sigma points. These sigma points are propagated through the true nonlinear system, without approximation, and a weighted mean and covariance are calculated. Before designing the estimator based on the UKF technique, first the extended Hovorka’s model 3 is discretized based on 5 min sampling time Ts Xk′+ 1 = Xk′ + Ts × f ′(Xk′ , uk) + Gk ωk , ωk ≈ N (0, Q ) yk = h(Xk′) + νk , νk ≈ N (0, R ) (6)

The details for the UKF algorithm with additive noise are given in the Table 3.30,36,37 For this estimator, the upper and lower bounds for states and the maximum possible changes of tmax,I and ke at each sampling 9849


Article


dead band is used as an alternative to CDEKF and UKF in this study. The main challenge with the MHE optimization approach is the potentially higher amount of computation effort required due to the complexity and size of the problem. Commercial and academic software packages have been developed to meet this challenge. APMonitor Modeling Language is used in this paper to implement the moving horizon estimation technique.43 In this method, the original Hovorka’s model is embedded in the MHE formulation and the tmax,I, ke, and UG parameters are considered to be time-varying parameters to be estimated every 5 min based on the recent measurement values obtained from the CGM sensor. The upper and lower bounds for states and maximum possible changes for time-varying parameters are considered to be the same as in UKF and CDEKF implementations (eq 5). The formulation for the MHE minimizing an L1-norm objective function with a dead band are

Table 3. Unscented Kalman Filtering equation χxk,i = [Xk̂ ′− 1

Xk̂ ′− 1 + Γ Pxk − 1

x χx− k,i = f(χk−1,i,uk−1) − 2n ̂ X′k = Σi = 0Wxi χx− k,i P−xk = Σi 2n= 0Wci (χx− k,i

constraints imposed

Xk̂ ′− 1 − Γ Pxk − 1 ]

̂−T − X̂ ′k−)(χx− k,i − X′k ) + Qk

γk,i = h(χx− k,i ) ŷk = Σi 2n= 0Wxi γk,i Pykyk = Σi2n= 0Wci (γk,i

yes yes yes yes yes

− ŷk)(γk,i − ŷk)T + Rk T ̂− Pxkyk = Σi2n= 0Wci (χx− k,i − X′k )(γk,i − ŷk) Kk = PxkykP−1 ykyk X̂ k′ = X̂ k′− + Kk(yk − ŷk)

yes

Pxk = P−xk − KkPykykKTk

α=1 β =2 κ=0 λ = α2(n + κ) − n

min Φ = WmT(eU + eL) + WpT(cU + cL) + ΔpT cΔp

X ,y ,p,d

⎛ dX ⎞ s . t . f ⎜ , X , y , p , d , u⎟ = 0 ⎝ dt ⎠

Γ= n+λ

Wxi =

λ n+λ λ + n+λ 1 2(n + λ)

Wci

1 2(n + λ)

Wx0 = Wc0 =

=

g (X , y , p , d , u) = 0

(1 − α 2 + β)

h(X , y , p , d , u) ≥ 0 db −y 2 db eU ≥ y − yx + 2 cL ≥ y ̂ − y eL ≥ yx −

time are identical to the CDEKF filter (eq 5). The constraints are imposed in different steps of the UKF algorithm as shown in Table 3.30 In the UKF algorithm, the constraints can be implemented more efficiently than the EKF algorithm by propagating the constraints through both the state and the covariance calculations. Moving Horizon Estimation. MHE is another widely used technique in dynamic systems to estimate state variables as well as constant or time-varying parameters. This technique calculates the states through the formulation of the estimation problem as an optimization problem incorporating the system model. MHE can also explicitly account for the physical constraints on the system by solving a constrained nonlinear optimization problem. MHE uses a window of prior measurements to minimize the deviation between model predictions and measurements by adjusting the parameters and states over a time horizon.38−44 Generally, Kalman filter-based techniques suffer from a number of limitations; for instance, physical constraints cannot be readily enforced, and poor quality data may result in the inaccurate estimation of the states. It is also difficult to obtain good estimates for the initial conditions of the error covariance matrix P, the covariance of the process noise Q, and the covariance of the observation noise R. For nonlinear or constrained systems, an optimization-based state estimation techniques such as the MHE approach are better suited for the estimation of the true system state. In this work, a MHE formulation is used that minimizes the L1-norm of the disturbance and output error variables. An important advantage of using an L1-norm objective function in comparison to the traditional least-squares form objective cost function in typical MHE formulations is the lower sensitivity of the L1-norm cost function to outliers, noise, or drifts in measurements.45−47The dead band in the objective is desirable for noise rejection, minimizing unnecessary parameter adjustments and movement of manipulated variables. Thus, the MHE technique minimizing the L1-norm objective function with a

cU ≥ y − y ̂ eL , eU , cL , cU ≥ 0

(7)

The nomenclature of the MHE algorithm is provided in Table 4. Table 4. Nomenclature for the MHE Minimizing an L1-Norm Objective Function with a Dead Band43 Φ yx y ŷ Wm Wp cΔp db x, u, p, d Δp f, g, h eL, eU cL , c U

objective function measurements (yx,0,..., yx,n)T model values (y0,..., yn)T prior model values (ŷ0,..., ŷn)T measurement deviation penalty penalty from the prior solution penalty from the prior parameter values dead band for noise rejection states (x), inputs (u), parameters (p), or unmeasured disturbances (d) change in parameters equation residuals, output function, and inequality constraints slack variable below and above the measurement dead-band slack variable below and above a previous model value

Partial Least Squares Models for PIC Estimator Individualization. The time-varying model parameters tmax,I and ke have a significant effect on the estimation of PIC values, and their appropriate initialization is crucial for improving the performance and convergence of the estimator. Because the values of these parameters differ from one patient to another, individualizing the PIC estimators to each patient is necessary. 9850


Article

Industrial & Engineering Chemistry Research The most accurate way to find these parameters is to fit the insulin compartmental model to the clinical PIC data for each patient. However, there are several limitations to this approach. For instance, accurate and reliable sensors that can report the PIC in real time do not exist, as the PIC values need to be measured by specific assays in clinical centers through an expensive and time-consuming operation. To the best of the authors’ knowledge, there is no existing method in the literature for the individualization of the PIC estimators. In this study, partial least-squares models are applied to obtain a good initial guess (a close approximation to the values obtained through model fitting) of these parameters. In reality, the demographic information on patients, such as weight, height, BMI, age, total daily insulin, and duration with T1DM, is readily available. Therefore, in this study we exploit this readily available information to determine a relationship between these demographic variables and the model parameters. Data-driven approaches, such as PLS regression models, are used when the precise mathematical relationship between two sets of data, demographic inputs (matrix V) and the output parameters (matrix Z), is not known. PLS is a multivariate regression method for modeling the relationship between two groups of data (explaining the variations in Z by using V) consisting of many noisy, correlated variables while handling potentially incomplete measurements with missing data.48 The mathematical relationship is derived from data collected from several experiments. The latent variables of the PLS model are identified to maximize the prediction performance of the model. This is achieved by finding components that maximize the covariance between the independent (V) and the dependent variables (Z). The V matrix consists of the demographic information on the patients, such as bodyweight, height, BMI, and total daily insulin dose, and the Z matrix is defined to be the two critical parameters tmax,I and ke that were determined for each patient using real PIC measurements. In this study, the estimators are evaluated based on two different initial values for the parameters tmax,I and ke, called training and testing parameters. The training parameters are obtained by fitting the insulin compartmental model to PIC data for each experiment. The testing parameters are computed in such a way that for each clamp study (with or without the IISWD) one experiment is selected and the PLS model is built based on the training parameters of the other experiments. Then the identified PLS model provided testing parameters which are tmax,I and ke for the excluded experiment. In addition to the parameters tmax,I and ke, the initial condition for UG should be defined. The value for UG is obtained from a steady-state assumption at initialization. Once the initial values for the parameters are determined from the PLS models developed using the training data, the proposed estimation techniques can be appropriately initialized for rapid convergence to provide good estimates of the PIC using only the CGM measurements and infused insulin data. Subjects and Clinical Study Experiments. Thirteen data sets from nine different subjects with T1DM are used.49 The study was conducted in adolescents with T1DM who attended the Yale Children’s Diabetes Clinic (New Haven, CT). Table 5 shows the demographic information on all nine subjects. Five-hour euglycemic clamps were conducted on two separate occasions for each subject separated by less than 8 weeks. One clamp was performed with the insulin infusion site warming device (IISWD) and the other without the IISWD. All subjects received insulin Aspart at 0.2 U/kg of body of weight bolus with

Table 5. Subject’s Demographic Information age (year) body weight (kg) height (cm) BMI (kg m−2) total daily insulin dose (U) (basal) duration of time with diabetes (month) HbA1c (%)

18.89 ± 1.62 60.56 ± 9.58 168.22 ± 9.19 21.39 ± 3.17 54.58 (24.71) ± 20.77 (11.87) 72.44 ± 44.88 7.36 ± 0.61

or without the IISWD at the start of the study. In all studies, the basal infusion of insulin via the insulin pump was suspended after the bolus was given, and a variable rate of 20% dextrose was infused and adjusted every 5 min based on bedside measurements of plasma glucose to maintain blood glucose levels within the range 90−100 mg/dL during the study.

■

RESULTS The performance of the proposed methods is tested using clinical data. The algorithms use only CGM and infused insulin data for estimation of PIC. The real PIC data are used only for comparison. In order to evaluate the performance of the proposed techniques, different performance indices have been used. The metrics root-mean-square error (RMSE) and mean absolute error (MAE) are calculated for a quantitative comparison of the estimation results The MAE measures the average magnitude of the absolute estimation errors over the test sample where all individual differences have equal weight. The RMSE also measures the average magnitude of the error, though using the square root of the average of squared differences between the estimated values and the actual data, which gives a relatively higher weight to the larger errors. Both MAE and RMSE provide the average estimation error in the same units as the variable of interest. The RMSE and MAE values are calculated as follows RMSE =

MAE =

ΣiN= 1(Ii − Iî)2 N

ΣiN= 1(|Ii − Ii|̂ ) N

(8)

(9)

where Iî is the estimated PIC and Ii is the measured PIC for the ith sample of N total measurement samples. The estimators are evaluated based on the training and testing parameters. In order to obtain the testing parameters for the clamp study with IISWD, a PLS model with four components was trained using weight, height, BMI, and duration with T1DM as inputs and the values of tmax,I parameters as the outputs. A PLS model with two components was trained with total daily insulin and duration with T1DM as inputs and values of ke parameters as outputs. For the clamp study without IISWD, the PLS model with four PLS components and tmax,I parameters values as outputs was built by weight, height, age, and duration with T1DM as inputs. The PLS models with four PLS components and ke as the output were developed by weight, height, BMI, and total daily insulin as inputs. The PLS model identified for predicting a particular parameter uses certain demographic information as inputs. These input variables are selected such that the average relative percentage error in predicting the testing parameters is minimized. Specifically, in each clamp study (with or without using IISWD) and for each parameter all available demographic information for each patient is considered initially. Then various 9851


Article


for the clamp with the IISWD and 9, 11.12, and 8.97 mU/L for the clamp without the IISWD, respectively. The values of these performance indices indicate significantly high correlations between the estimated and the measured insulin concentrations based on training parameters. Tables S3 and S4 show the RMSE and MAE values for testing parameters of both clamp studies (with and without the IISWD) based on CDEKF, UKF, and MHE. The average RMSE values are found to be 26.90, 29.03, and 27.02 mU/L based on CDEKF, UKF, and MHE for the clamp with the IISWD and 19, 13.14, and 18.04 mU/L for the clamp without the IISWD, correspondingly. The average MAE values are found to be 23.04, 24.56, and 22.97 mU/L based on CDEKF, UKF, and MHE for the clamp with the IISWD and 14.83, 10.56, and 14.24 mU/L for the clamp without the IISWD, respectively. The values of these performance indices indicate significantly high correlations between the estimated and the measured insulin concentrations based on testing parameters as well. Figure 3 compares the estimated PIC values with the real measurements for both training and testing parameters of experiment 3 with the IISWD based on UKF, CDEKF, and MHE. Figure 4 compares the estimated PIC values with the real measurements for both training and testing parameters of experiment 5 without the IISWD based on UKF, CDEKF, and MHE. Figures 5 and 6 compare the estimated CGM with the real CGM measurements for experiments 3 with the IISWD and 5 without IISWD based on testing and training parameters to illustrate the satisfactory performance of CDEKF, UKF, and MHE estimators in tracking and estimation of model output. The proposed estimators can track and estimate the CGM data well and can accurately estimate the PIC.

PLS models are identified using different subsets of the demographic variables as inputs. Finally, the subset of the demographic variables that provides the best prediction of the parameters (minimum prediction error) is reported as the most accurate PLS models (Table 6). Hence, in each clamp study a Table 6. PLS Models for PIC Estimator Individualization clamp study

tmax,I

with IISWD weight, height, BMI, and duration with T1DM without weight, height, age, and duration IISWD with T1DM

ke total daily insulin and duration with T1DM weight, height, BMI, and total daily insulin

different subset of the demographic variables provided the best prediction results, which is expected since the pharmacokinetics of insulin change when using the infusion site warming device. The average relative error percentages of testing parameters for ke parameter in the clamp study with and without IISWD are 17.22% and 13.66%, respectively. For the tmax,I parameter in the clamp study with and without IISWD, the average relative errors of testing parameters are found to be 2.35% and 3.46%, respectively. Figures 1 and 2 compare the estimated PIC values with the real measurements for seven experiments under the clamp study with the IISWD and six experiments under the clamp study without the IISWD using CDEKF, respectively. Tables S1 and S2 show the RMSE and MAE values for training parameters of both clamp studies (with and without the IISWD) based on CDEKF, UKF, and MHE. The average RMSE values are found to be 13.22, 13.38, and 14.65 mU/L based on CDEKF, UKF, and MHE for the clamp with the IISWD and 12.13, 13.83, and 12.10 mU/L for the clamp without the IISWD, correspondingly. The average MAE values are found to be 9.84, 10.18, and 11.02 mU/L based on CDEKF, UKF, and MHE

Figure 1. Comparison of estimated and measured PIC for the clamp study with the IISWD based on CDEKF: (blue line) estimated PIC based on training parameters, (red line) estimated PIC based on testing parameters, (filled green circles) measured PIC. 9852


Article


Figure 2. Comparison of estimated and measured PIC for the clamp study without the IISWD based on CDEKF: (blue line) estimated PIC based on training parameters, (red line) estimated PIC based on testing parameters, (filled green circles) measured PIC.

Figure 3. Comparison of estimated and measured PIC data based on CDEKF, UKF, and MHE for experiment 3 with the IISWD: (blue line) estimated PIC based on training parameters, (red line) estimated PIC based on testing parameters, (filled green circles) measured PIC.

■

experiments are higher because of these fluctuations. The fluctuations might be due to the measurement errors of the assay used for PIC measurement. Another possible reason for these abrupt changes might be due to subcutaneous insulin kinetics. In this case, Hovorka’s model cannot describe these fluctuations well, and consequently, the estimators have relatively lower

DISCUSSION

The performance of the proposed methods is relatively poorer for two experiments as is readily observed in the above tables (experiments six and four with and without the IISWD, respectively). There are high fluctuations in the measured PIC values in these experiments. The RMSE and MAE values for both 9853


Article


Figure 4. Comparison of estimated and measured PIC data based on CDEKF, UKF, and MHE for experiment 5 without the IISWD: (blue line) estimated PIC based on training parameters, (red line) estimated PIC based on testing parameters, (filled green circles) measured PIC.

Figure 5. Comparison of estimated and measured CGM data based on CDEKF, UKF, and MHE for experiment 3 with the IISWD: (blue line) estimated CGM based on training parameters, (red line) estimated CGM based on testing parameters, (filled green circles) measured CGM.

estimation accuracy for such undefined fluctuations. Despite the high performance of estimator for training parameters of experiment 7 with the IISWD, the poor PIC estimation was observed for testing parameters. For this case, the PLS model built could not provide accurate testing parameters. This performance may be due to the significant variance in duration of time with diabetes, which ranges from 24 to 144 months. This

variability in a limited number of subjects (6 for clamp study with IISWD) can cause inaccuracies in the model predications. Nevertheless, the model prediction results demonstrate poor predictions for one experiment, which may be attributed to the insufficient data available for training the PLS models. Identifying the PLS models with a greater number of subjects and 9854


Article


Figure 6. Comparison of estimated and measured CGM data based on CDEKF, UKF, and MHE for experiment 5 without the IISWD: (blue line) estimated CGM based on training parameters, (red line) estimated CGM based on testing parameters, (filled green circles) measured CGM.

allows the analysis of the performance under intrasubject variability, including physical activity. Incorporating the effect of physical activity in plasma insulin concentration dynamics would lead to improved performance of an artificial pancreas control system during exercise. As the plasma insulin concentrations measurements are not available during exercise in the clinical data sets used in this work, the IISWD-based results can be considered as a surrogate. The effects of the insulin infusion site warming device on patients was analyzed for individuals included in both with and without IISWD experiments. Using the IISWD accelerates the insulin action through increased blood flow to the area of insulin absorption and correspondingly the tmax,I (time-to-maximum of the absorption of subcutaneously injected insulin) decreases. On the other hand, there is no consistent change in the parameter ke between the experiments with and without infusion site warming device observed in the results. The unique properties of the proposed techniques which make them superior to the other methods for determining the active insulin in the body are their ability to first be individualized by using PLS models built using the demographic information on the subjects and then capture the changes in glucose−insulin dynamics by estimating effective model parameters over time. In this study, two effective parameters tmax,I and ke on PIC and one parameter UG for gut absorption rate are considered to be time varying in the original model to describe the temporal system dynamics. To show the importance of individualization of PIC estimators and estimation of effective model parameters, the average RMSE and MAE of estimated PIC values were also calculated for the estimators designed based on the original Hovorka’s model. The average RMSE were found to be 38.17 and 46.73 based on the CDEKF, 38.39 and 50.67 based on UKF, and 40.17 and 49.84 based on MHE for the data sets with and without IISWD, respectively. The average MAE were found to be 33.06

encompassing a wide range of demographic information would significantly improve the accuracy of the predication. A requirement for identifying highly accurate PLS models is the availability of rich data with an abundant number of training samples from numerous clinical experiments. To train the PLS model per patient to model intrasubject variability, a sufficient quantity of data for each patient is necessary. The clinical experiments considered in this work consist of two data sets for some of the subjects. Building data-driven models, such as PLS models, per patient to characterize the intrasubject variability is the goal of future studies where sufficient quantity of data are collected for each patient. There are different metrics to consider the effect of active insulin in an AP control system. In this work, one of these metrics, the PIC, is investigated. Providing PIC information can reduce hypoglycemia by preventing the AP system from administering additional insulin when the PIC is already high. However, insulin is measured only from intravenous blood samples and analyzed off-line, which is not suitable for an AP system that requires this information in real time. The proposed PIC estimators based on CDEKF, UKF, and MHE methods can provide an estimation of PIC in real time. These methods use only CGM measurements and infused insulin data in real time and do not require any manual inputs such as meal information from the patients. The results show that the proposed methods have satisfactory performance in estimating PIC. In this study, PIC after only a single bolus insulin infusion is estimated. The proposed methods are now incorporated into the integrated multivariable adaptive AP system, and clinical trials are being conducted to test these estimators over longer periods of time with multiple insulin boluses and basal infusions.3,8 Exercise and the use of an infusion site warming device (IISWD) both have a similar effect on insulin pharmacokinetics since both mechanisms increase the body temperature, which 9855



■

and 39.47 based on the CDEKF, 33.10 and 42.66 based on UKF, and 34.38 and 42.67 based on MHE for the data sets with and without IISWD, respectively. All these average RMSE and MAE values prove the significant role of an adaptive individualized PIC estimator. However, variations in insulin−glucose dynamics because of differences in daily activities or metabolic state of the patients may affect additional parameters of the model. The proposed methods can also take into account the effects of these factors by including more model parameters as extended states (in CDEKF and UKF) or as manipulated variables (in MHE) to the original model. The observability of the augmented system with the addition of new parameters as extended states should be confirmed. If the augmented system with additional parameters as extended states is observable, an observer can be designed for the new augmented system to simultaneously estimate both the states and the parameters. This is a necessary condition for the real-time estimation of the time-varying model parameters. An important consideration in the design of AP systems is the online computational tractability of the control and estimation algorithms. For the proposed PIC estimators, CDEKF has the lowest computational requirements, while UKF and MHE have greater computational complexity. The MHE approach requires the greatest computational effort since a nonlinear programming problem must be solved at each sampling instance. However, decreasing the computational complexity may adversely affect the prediction accuracy as these are typically competing criteria. For online and real-time applications, it is necessary to consider the accuracy of the estimated PIC values along with the associated computational burden, especially when implementing the algorithms on computationally constrained wearable or implantable hardware. In this study, we design and adapt efficient and widely used observers for estimation of PIC. Nevertheless, each of the applied estimation techniques have both advantages and disadvantages in practice. Our end goal is to test these techniques in clinical experiments under different situations and evaluate their performance.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: (312) 567-3042. Fax: (312) 5678874. ORCID

Ali Cinar: 0000-0002-1607-9943 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS

■

REFERENCES

This work was supported by the National Institutes of Health (NIH) under grants 1DP3DK101075-01 and 1DP3DK10107701.

(1) El-Khatib, F. H.; Russell, S. J.; Nathan, D. M.; Sutherlin, R. G.; Damiano, E. R. Bi-hormonal closed-loop blood glucose control for type 1 diabetes. Sci. Transl. Med. 2010, 2, 27ra27. (2) Steil, G. M.; Rebrin, K.; Darwin, C.; Hariri, F.; Saad, M. F. Feasibility of automating insulin delivery for the treatment of type 1 diabetes. Diabetes 2006, 55, 3344−3350. (3) Turksoy, K.; Quinn, L. T.; Littlejohn, E.; Cinar, A. An integrated multivariable artificial pancreas control system. J. Diabetes Sci. Technol. 2014, 8, 498−507. (4) Boiroux, D.; Duun-Henriksen, A. K.; Schmidt, S.; Nørgaard, K.; Poulsen, N. K.; Madsen, H.; Jørgensen, J. B. Adaptive control in an artificial pancreas for people with type 1 diabetes. Control Eng. Pract. 2017, 58, 332−342. (5) Cameron, F.; Ly, T.; Forlenza, G.; Patek, S.; Baysal, N.; Messer, L.; Clinton, P.; Maahs, D.; Buckingham, B.; Bequette, B. Inpatient clinical trial of a fully closed-loop artificial pancreas using only CGM and accelerometer data for insulin dosing. Diabetes Technol. Ther. 2016, A23−A24. (6) Kovatchev, B.; Tamborlane, W. V.; Cefalu, W. T.; Cobelli, C. The artificial pancreas in 2016: a digital treatment ecosystem for diabetes. Diabetes Care 2016, 39, 1123−1126. (7) Trevitt, S.; Simpson, S.; Wood, A. Artificial pancreas device systems for the closed-Loop control of type 1 diabetes what systems are in development? J. Diabetes Sci. Technol. 2016, 10, 714−723. (8) Cinar, A.; Turksoy, K.; Hajizadeh, I. Multivariable artificial pancreas method and system. U.S. Patent Appl. US15/171,355, 2016. (9) Nixon, R.; Pickup, J. C. Fear of hypoglycemia in type 1 diabetes managed by continuous subcutaneous insulin infusion: is it associated with poor glycemic control? Diabetes Technol. Ther. 2011, 13, 93−98. (10) Ellingsen, C.; Dassau, E.; Zisser, H.; Grosman, B.; Percival, M. W.; Jovanovic, L.; Doyle, F. J. Safety constraints in an artificial pancreatic β cell: an implementation of model predictive control with insulin on board. J. Diabetes Sci. Technol. 2009, 3, 536−544. (11) Toffanin, C.; Zisser, H.; Doyle, F. J.; Dassau, E. Dynamic insulin on board: incorporation of circadian insulin sensitivity variation. J. Diabetes Sci. Technol. 2013, 7, 928−940. (12) Steil, G. M.; Palerm, C. C.; Kurtz, N.; Voskanyan, G.; Roy, A.; Paz, S.; Kandeel, F. R. The Effect of insulin feedback on closed loop glucose control. J. Clin. Endocrinol. Metab. 2011, 96, 1402−1408. (13) Revert, A.; Garelli, F.; Picó, J.; De Battista, H.; Rossetti, P.; Vehí, J.; Bondia, J. Safety auxiliary feedback element for the artificial pancreas in type 1 diabetes. IEEE Trans. Biomed. Eng. 2013, 60, 2113−2122. (14) Palerm, C. C. Physiologic insulin delivery with insulin feedback: a control systems perspective. Comput. Methods Programs Biomed. 2011, 102, 130−137. (15) Ruiz, J. L.; Sherr, J. L.; Cengiz, E.; Carria, L.; Roy, A.; Voskanyan, G.; Tamborlane, W. V.; Weinzimer, S. A. Effect of insulin feedback on closed-loop glucose control: a crossover study. J. Diabetes Sci. Technol. 2012, 6, 1123−1130.

■

CONCLUSION Incorporating information about the insulin-on-board is crucial for AP systems to calculate the optimum insulin infusion rate. Knowing the PIC values can provide valuable information on IOB. PIC cannot be measured in real time for use in AP systems, though it can be accurately inferred through estimation algorithms based on reliable and observable glucose−insulin dynamic models. In this study, the CDEKF, UKF, and MHE approaches are proposed to provide an accurate estimate of PIC based on CGM data collected every 5 min and infused insulin data. The proposed methods are able to estimate PIC in real time satisfactorily. These methods would be beneficial in an AP system for the real-time estimation of unmeasurable variables such as plasma insulin concentrations.

■

Article

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.iecr.7b01618. Performance evaluation of the proposed algorithms for the clamp study with the IISWD based on training parameters, without the IISWD based on training parameters, with the IISWD based on testing parameters, and without the IISWD based on testing parameters (PDF) 9856


Article

Industrial & Engineering Chemistry Research (16) León-Vargas, F.; Garelli, F.; De Battista, H.; Vehí, J. Postprandial response improvement via safety layer in closed-loop blood glucose controllers. Biomed. Signal Process. Control 2015, 16, 80−87. (17) de Pereda, D.; Romero-Vivo, S.; Ricarte, B.; Rossetti, P.; Ampudia-Blasco, F. J.; Bondia, J. Real-time estimation of plasma insulin concentration from continuous glucose monitor measurements. Comput. Methods Biomech. Biomed. Eng. 2016, 19, 934−942. (18) Eberle, C.; Ament, C. The unscented Kalman Filter estimates the plasma insulin from glucose measurement. BioSystems 2011, 103, 67− 72. (19) Eberle, C.; Ament, C. Real-time state estimation and long-term model adaptation: a two-sided approach toward personalized diagnosis of glucose and insulin levels. J. Diabetes Sci. Technol. 2012, 6, 1148− 1158. (20) Kiriakidis, K.; O’Brien, R. Optimal estimation of blood insulin from blood glucose. International Mechanical Engineering Congress and Exposition; ASME, 2006; pp 1175−1180. (21) Kiriakidis, K. Particle filter for plasma insulin estimation. Dynamic Systems and Control Conference; ASME, 2010; pp 369−372. (22) Kiriakidis, K.; O’Brien, R. Estimation of plasma insulin using nonlinear filtering. International Mechanical Engineering Congress and Exposition; ASME, 2007; pp 203−208. (23) Hajizadeh, I.; Turksoy, K.; Cengiz, E.; Cinar, A. The extended and unscented Kalman filtering methods for real-time plasma insulin concentration estimation in an artificial pancreas control system for patients with type 1 diabetes. Annual Meeting AIChE; AIChE, 2016. (24) Hajizadeh, I.; Turksoy, K.; Cengiz, E.; Cinar, A. Real-time estimation of plasma insulin concentration using continuous subcutaneous glucose measurements in people with type 1 diabetes. Proceedings of the American Control Conference (ACC), Seattle, WA; ACC, 2017. (25) Hovorka, R.; Canonico, V.; Chassin, L. J.; Haueter, U.; MassiBenedetti, M.; Federici, M. O.; Pieber, T. R.; Schaller, H. C.; Schaupp, L.; Vering, T.; Wilinska, M. E. Nonlinear model predictive control of glucose concentration in subjects with type 1 diabetes. Physiol. Meas. 2004, 25, 905. (26) Dalla Man, C.; Rizza, R. A.; Cobelli, C. Meal simulation model of the glucose-insulin system. IEEE Trans. Biomed. Eng. 2007, 54, 1740− 1749. (27) Bergman, R. N.; Ider, Y. Z.; Bowden, C. R.; Cobelli, C. Quantitative estimation of insulin sensitivity. Am. J. Physiol. Endocrinol. Metab. 1979, 236, E667. (28) Dochain, D. State and parameter estimation in chemical and biochemical processes: a tutorial. J. Process Control 2003, 13, 801−818. (29) Crassidis, J. L.; Junkins, J. L. Optimal estimation of dynamic systems; CRC Press, 2011. (30) Kolås, S.; Foss, B.; Schei, T. Constrained nonlinear state estimation based on the UKF approach. Comput. Chem. Eng. 2009, 33, 1386−1401. (31) Haseltine, E. L.; Rawlings, J. B. Critical evaluation of extended Kalman filtering and moving-horizon estimation. Ind. Eng. Chem. Res. 2005, 44, 2451−2460. (32) Simon, D. Optimal state estimation: Kalman, H infinity, and nonlinear approaches; John Wiley & Sons, 2006. (33) Simon, D.; Simon, D. L. Kalman filtering with inequality constraints for turbofan engine health estimation. IEE Proc.: Control Theory Appl. 2006, 153, 371−378. (34) Kandepu, R.; Foss, B.; Imsland, L. Applying the unscented Kalman filter for nonlinear state estimation. J. Process Control 2008, 18, 753−768. (35) Simon, D.; Simon, D. L. Aircraft turbofan engine health estimation using constrained Kalman filtering. J. Eng. Gas Turbines Power 2005, 127, 323−328. (36) Julier, S. J.; Uhlmann, J. K.; Durrant-Whyte, H. F. A new approach for filtering nonlinear systems. Proceedings of the American Control Conference; ACC, 1995; pp 1628−1632. (37) Van Der Merwe, R. Sigma-point Kalman filters for probabilistic inference in dynamic state-space models. Ph.D. thesis, Oregon Health & Science University, 2004.

(38) Hashemian, N.; Armaou, A. Simulation, model-reduction, and state estimation of a two-component coagulation process. AIChE J. 2016, 62, 1557−1567. (39) Hashemian, N.; Armaou, A. Fast moving horizon estimation of nonlinear processes via Carleman linearization. American Control Conference (ACC); ACC, 2015; pp 3379−3385. (40) Safdarnejad, S. M.; Hedengren, J. D.; Lewis, N. R.; Haseltine, E. L. Initialization strategies for optimization of dynamic systems. Comput. Chem. Eng. 2015, 78, 39−50. (41) Safdarnejad, S. M.; Hedengren, J. D.; Baxter, L. L. Dynamic optimization of a hybrid system of energy-storing cryogenic carbon capture and a baseline power generation unit. Appl. Energy 2016, 172, 66−79. (42) Safdarnejad, S. M.; Gallacher, J. R.; Hedengren, J. D. Dynamic parameter estimation and optimization for batch distillation. Comput. Chem. Eng. 2016, 86, 18−32. (43) Hedengren, J. D.; Shishavan, R. A.; Powell, K. M.; Edgar, T. F. Nonlinear modeling, estimation and predictive control in APMonitor. Comput. Chem. Eng. 2014, 70, 133−148. (44) Eaton, A. N.; Safdarnejad, S. M.; Hedengren, J. D.; Moffat, K.; Hubbell, C. B.; Brower, D. V.; Brower, A. D. Post-installed fiber optic pressure sensors on subsea production risers for severe slugging control. 34th International Conference on Ocean, Offshore and Arctic Engineering, Newfoundland, Canada, 2015; ASME, 2015; p V05BT04A055. (45) Safdarnejad, S. M.; Hedengren, J. D.; Baxter, L. L. Plant-level dynamic optimization of cryogenic carbon capture with conventional and renewable power sources. Appl. Energy 2015, 149, 354−366. (46) Safdarnejad, S. M.; Hedengren, J. D.; Baxter, L. L.; Kennington, L. Investigating the impact of cryogenic carbon capture on power plant performance. Proceedings of the American Control Conference (ACC), Chicago, IL; ACC, 2015; pp 5016−5021. (47) Hedengren, J. D.; Eaton, A. N. Overview of estimation methods for industrial dynamic systems. Optim. and Eng. 2017, 18, 155−178. (48) Cinar, A.; Palazoglu, A.; Kayihan, F. Chemical process performance evaluation; CRC Press, 2007. (49) Cengiz, E.; Weinzimer, S. A.; Sherr, J. L.; Tichy, E. M.; Carria, L.; Cappiello, D.; Steffen, A.; Tamborlane, W. V. Faster in and faster out: accelerating insulin absorption and action by insulin infusion site warming. Diabetes Technol. Ther. 2014, 16, 20−25.

9857


Plasma Insulin Estimation in People with Type 1 Diabetes Mellitus

Recommend Documents