Online Frequency Domain Volterra Model of Glucose-Insulin Process in Type-1 Diabetics
1Department of Electrical Engineering, Bengal Engineering and Science University, Shibpur, Howrah, India
Modern close loop control for blood glucose level in a diabetic patient necessarily uses an explicit model of the process. A fixed parameter full order or reduced order model does not characterize the inter-patient and intra-patient parameter variability. This paper deals with an online frequency domain kernel estimation method for modeling a nonlinear dynamic system of multivariable glucose-insulin process in a type-1 diabetic patient that captures the process dynamics in presence of uncertainties and parameter variations. The present work proposes a frequency domain kernel estimation of a Volterra model using the harmonic excitation input by taking FFT on the input data sequence from the glucose-insulin process of the patient. Volterra equations up to second order kernels with extended input vector for Volterra model are solved online by adaptive recursive least square (ARLS) algorithm. Twice the length of the extended input vector for the glucose-insulin process is considered for finding the frequency domain kernels that can be directly used as the Volterra transfer function and are useful for closed loop internal model control. The input-output data taken from the 19th order first principle model of the patient in intravenous route, have been used to identify the system with a short filter memory length of M=2 and the validation results have shown good fit both in frequency and time domain responses with nominal patient as well as with intrapatient parameter variations.
At a glance: Figures
Keywords: System identification, nonparametric model, glucose-insulin interaction, Volterra model, frequency domain kernels, Volterra transfer function
Journal of Biomedical Engineering and Technology, 2014 2 (2),
Received February 28, 2014; Revised May 09, 2014; Accepted May 12, 2014Copyright © 2013 Science and Education Publishing. All Rights Reserved.
Cite this article:
- Bhattacharjee, A., and A. Sutradhar. "Online Frequency Domain Volterra Model of Glucose-Insulin Process in Type-1 Diabetics." Journal of Biomedical Engineering and Technology 2.2 (2014): 13-20.
- Bhattacharjee, A. , & Sutradhar, A. (2014). Online Frequency Domain Volterra Model of Glucose-Insulin Process in Type-1 Diabetics. Journal of Biomedical Engineering and Technology, 2(2), 13-20.
- Bhattacharjee, A., and A. Sutradhar. "Online Frequency Domain Volterra Model of Glucose-Insulin Process in Type-1 Diabetics." Journal of Biomedical Engineering and Technology 2, no. 2 (2014): 13-20.
|Import into BibTeX||Import into EndNote||Import into RefMan||Import into RefWorks|
In fast lifestyle pattern of these days, intensive treatment of type-1 diabetic patients is required to avoid later life complications. Proper dose of insulin is applied to the bloodstream intermittently through subcutaneous or intravenous infusion to regulate the blood glucose (BG) level. A tight closed loop control of BG is required to maintain normoglycaemia (BG level in the range of 63-144 mg/dl) in presence of normal meal and activity conditions of the patient [1-6,19]. Effective control of biological processes is a non-trivial task requiring a model that is able to adequately acquire the dynamic behaviour of the process over its complete operating range. The complex nonlinear process of glucose metabolism is linked to a number of internal factors. With occasional blood glucose sensing, routine food intake and other activity conditions, the system appears highly stochastic. The closed loop control involves interplay between the nonlinear dynamics of the physiological process, the inter- patient and intra-patient variability, uncertainties and disturbances in implantable infusion device (actuator) and the glucose sensor [7, 8]. Researchers worked on various close loop control algorithms for BG regulation either based on first principle full order or Bergman’s minimal model of the glucose-insulin (GI) process [9-18]. Modern control methodologies though demonstrated adequate performance; the inherent uncertainty in the model (or patient) has not been explicitly addressed. This omission can lead to significant performance degradation if the effect of parameter variations is not reflected in the derived model. To avoid retuning of the controller every time for different patient conditions, the model should be capable of online capturing of such dynamics .
Though the theory of data driven identification techniques are well-established [19, 20] in last three decades, very limited work on modeling of Glucose-Insulin (GI) process in frequency domain has been reported. On-line ARX model and nonparametric time domain kernels of the GI process have been reported elsewhere by the present authors [21, 22]. Nonlinear ARX (NARX) model is significant among different black box model structures used to capture nonlinear dynamics but determination of order and structure of such model on-line is a difficult task even for a single input single output (SISO) system and the problem is further compounded for multi-variable systems . The parametric models also lack full information of nonlinearities in the coupled systems. Nonparametric model, on the contrary, decomposes arbitrary basis functions leading to most concise signal representation and best conditions of identification algorithms [23, 24]. In a previous work by the present authors  the frequency domain kernel identification method was established and validated using Volterra- Hammerstein structure for the Guyton model of the patient. A potential advantage of using nonparametric model is that it can yield nonlinear model predictive control directly from the identified process using the Volterra transfer function of the process derived from the frequency domain kernels [25, 26, 27, 28].
In the present paper, an online nonparametric identification method has been presented where the frequency domain kernels are computed for a simpler Volterra-only structure that captures the dynamics of the GI process with parameter variations in type-1 diabetic patients. The 19th order first principle Sorensen model of the type-1 diabetic patient embedded into a hardware platform is used as the virtual patient to be identified. The important parameter of the patient is varied to realize different patient conditions. The resulting model has been validated for a wide range of input, meal disturbance and patient parameter variations.
2. System Overview
The closed loop drug delivery system for BG regulation necessarily requires a patient model estimator, which continually estimates patient states from the physiological process data in presence of disturbances and parameter variations from the input-output Data. The block diagram of the closed loop system is shown in Figure 1.
A pharmacokinetic and pharmacodynamic compartmental model of glucose-insulin (GI) process involving various organs has been constructed as given by Guyton et al. , that was modified by Sorensen  resulting in large number of system states to accommodate inter organ transmissions. Utilizing compartmental modeling techniques, the glucose-insulin model of diabetic p tient is represented schematically in Figure 2.
The body has been divided into six physiologic compartments as represented in  and  shown in Figure 2 where Brain- represents the central nervous system; Heart/Lung - represents the rapidly mixing vascular volumes of the heart, lungs and arteries; Periphery - includes skeletal muscle and adipose tissues; Gut -includes the effects of stomach and intestine; Liver – has two way function on glucose and Kidney for excretion. Arrows represents the direction of flow of blood in the various compartments. Mass balance equations around tissues important to glucose and insulin dynamics have been used. It has been assumed that the glucose or insulin concentration in a compartment is in equilibrium with the blood leaving the given compartment.
In the present work, a meal disturbance model is included in the model description as developed by Lehmann et al.  and Parker et al.  and includes uncertainty modeling [3, 18] (details given in Appendix). To maintain relatively constant carbohydrate release from the stomach during intestinal adsorption gastric emptying (Gempt) of carbohydrate is modeled. The rate of gastric emptying (Gempt), increases linearly from the instant of taking meal, then plateaus to a maximum value for some time; and ultimately decreases to zero linearly. The duration for which the gastric emptying rate remains constant at a maximum value is a function of the amount of carbohydrate intake. The rise to and fall from this maximum rate is a linear function (ramp), taking place over a thirty minute span. Small meals (<10.2 g carbohydrate).only contained the rise and fall phases (triangle) shape, never reaching the plateau emptying rate. Larger meals are described by a trapezoidal wave. A firstorder filter is used after this Gempt function for smoothening the output .
In the present work, intravenous (I/V) insulin delivery system is considered due to its inherent advantages of rapid delivery with negligible dead-time, a higher percentage of drug utilization, less amount of drug required, faster response to insulin over-delivery and potential for improved closed-loop controller performance.
3. Identification Algorithm of Nonlinear Model Structure
Successful identification requires a proper experimental design and choice of a proper model structure. Experimental design involves selection of experimental parameters like sampling time and excitation signal. Many researchers have used block oriented models that give insight to the complexity of the process, particularly for nonlinear processes [22-27]. Use of block structure gives an advantage of approximating the nonlinear processes with higher accuracy and enhances the chances to get a reduced order model, which is uncomplicated and computationally efficient. Nonparametric approaches represent the dynamic characteristics of the nonlinear system from expansion of the Volterra series or kernel.3.1. Volterra Model
The finite Volterra series up to second order kernel for a multivariable (MISO) system with xi inputs, where i=1, 2, … m the output y(t) for a kernel memory length M is expressed as:
The self-kernels hii acting on a single input are symmetric and the cross-kernels hij acting on different inputs and are asymmetric. The second order Volterra structure for the present single-input single output GI process thus follows directly from equation 2 and the output of the process is given by:
where, is the zero order Volterra kernel and , are the first and second order Volterra kernels respectively associated with the input x(t) in time domain and ’*’ denotes the convolutions. The output y(t) can also be expressed in terms of the nonlinear kernel operator H(t) as:
Since the higher order model is not very useful for computation and simulation purpose, the second order Volterra model having a short memory of M = 2 has been considered here. For the present system, Volterra equations up to second order kernels with extended input vector are solved online by adaptive recursive least square (ARLS) algorithm .
Adaptive recursive least squares (ARLS) algorithm has been used to find the filter coefficients recursively on the basis of least squares of the error signal. The structural diagram of the RLS filter is shown in Figure 3.
The idea behind RLS filters is to select the filter coefficients and update the same with new data set by minimizing the following cost function involving the error signal , desired output and model output :
Figure 3 shows the ARLS structure with instantaneous value of the variables. is a factor that controls the memory span of the adaptive filter. The filter coefficients are adapted according to the following steps:
I. Initialization: Define filter memory where is a small positive constant;
II. Operations: for t = 1 to number of iterations.
a) Create the extended input vector:
b) Compute the output:
c) Compute the error:
d)Compute the gain matrix:
e) Update the filter vector:
f) Update the inverse autocorrelation matrix P:
In the above representations, the functions represent the kernels associated to the nonlinear operators The above representations have the same memory for nonlinearity orders. In this case the second order Volterra kernel is a matrix. If we consider symmetric kernels of memory the second order Volterra kernel requires the determination of coefficients for single input. Hence the total number of coefficients is .3.3. Frequency Domain Representation of Volterra Model
Considerable amount of effort has been made in recent years to develop frequency domain representation of nonlinear systems based on Volterra series representation of input output behavior [25, 26, 27]. An online frequency domain kernel estimation method for Volterra model has been presented here for the dynamic system of the GI process. The main advantage of using frequency domain kernel is to get the Volterra Transfer Function (VTF) that can be directly used in model predictive control.
The input-output relationship in frequency domain can be easily obtained as follows. Let and denote the Fourier transforms of the inputs and output respectively. The Fourier transform of the order kernel (nonlinear impulse response of order k) can be written as
is called the nonlinear transfer function or Volterra Transfer Function (VTF) of order k, where is the maximum order of the kernel, which can be determined experimentally by applying a multi-tone signal at the input to the system and measuring the response. The output for the two input single output system is given by:
where represents the corresponding frequency shifts. is computed by using Overlap-Save algorithm. The frequency domain kernels have been computed by taking the FFTs on respective time domain kernels for a specific length of extended input vector.
The overlap-save algorithm reconstructs time domain output y(t) by taking inverse FFT of the output :
By the overlap-save algorithm, the convolution and correlation operations are very easily and efficiently implemented in frequency domain, which is rather computationally intensive when performed directly by the time-domain ARLS method.
The steps of the Overlap-Save algorithm are:
I. Let N be number of data points in time domain extended input vector. N zeros are added to the left of the input vector that contains both the insulin and meal inputs.
II. The length of the input data points is twice the length of the time domain kernels i.e. 2N. Time domain Volterra kernels are computed from 2N extended input vector, so that the result of the FFT will be the same length as that of the FFTs of the input sections.
III. 2N point FFT taken for both the input vector and the kernels and stored in memory.
IV. The two FFTs are now multiplied to generate the output in frequency domain (each element in one of the arrays is multiplied by the corresponding element in the other). This corresponds to convolution in time domain.
V. The first row of the result gives the output. The first half is added to an array as the output of the filter for the given input block.
VI. Finally the input block is updated.3.4. Flow Chart of the Identification Algorithm
4. Experiments and Verification of Results
As the physiological process in type-1 diabetic patient is completely lacking endogenous insulin secretion, the BG level rises to a large value exceeding 350mg/dl in open loop when 50gm meal is given orally. The use of pseudo random binary sequence (PRBS) of 10% over basal dose of insulin is a normal practice for offline identification of a process from its input-output data  but not practicable for online identification. Here for online identification of the data driven model of the GI process of a type-1 diabetic patient, a basal insulin dose of 22.33mU/min has been used. The patient is subjected to a prescribed 50gm meal (equivalent carbohydrate taken orally (CHO)) ingestion at 300min [2, 3]. The open loop response of blood glucose level with constant insulin infusion and meal disturbance for 900min span is shown in Figure 4.
The actual output and the identified model output generated by simulation for above input-output condition is shown in Figure 5. The response shows a good match between the actual glucose output and the output from the identified model for nominal patient parameter values.
Uncertainty may be related to variations in patient parameters. The Sorensen full order model  identifies the metabolic terms as most responsible for changes in glucose and insulin dynamics, and can be analyzed by the following threshold functions of glucose metabolism [3, 18]:
Here the subscript i is the state vector element involved in the metabolic effect and e denotes specific effects within the model, such as the effect of glucose on hepatic glucose uptake (EGHGU), or the effect of insulin on peripheral glucose uptake (EIPGU). Inter- or intra-patient variability is classified physiologically as either a receptor () parameter or post-receptor () parameter defect. This uncertainty formulation implies a structured effect of variability on the model, such that the tissues most important to parametric uncertainty are the liver and peripheral tissues . Differences in insulin clearance between patients also exist, and can be modeled as deviations in the fraction of hepatic insulin clearance (FHIC) [3, 18]. The parameters EGHGU- and EIPGU- are varied by ±40% from their nominal values of -5.82 and -1.48 respectively and FHIC is varied by ± 20% from its nominal value of 0.4 . The effect of intrapatient variability in FHIC in between meals is shown in Figure 6.
The validation was also tested using the power spectral densities (PSD) of the frequency domain output of the nominal patient model and that of the desired output of the system as plotted in Figure 7. This also shows very good matching between the frequencies where the power of the two signals lies. The closeness of these two PSDs also confirms the model’s accuracy.
Model-based predictive control of insulin delivery system requires development of an accurate model of the human glucose-insulin system and design of a constrained controller. A data driven block-oriented modeling in time domain as well as in frequency domain for the nonlinear dynamic system of multivariable glucose-insulin process in a type-1 patient with relatively short memory effects has been developed. The advantages of block-oriented model have been utilized with proper selection of Volterra kernels by ARLS algorithm and extended input vectors for the nonlinear process containing both deterministic insulin input as well as meal disturbances. Each time domain Volterra term is expressed by means of a multifold convolution integral and each frequency domain Volterra term is expressed by means of a product with FFT inputs over the finite memory interval. Good agreement was found between the response from the actual process and the identified model. The set of kernels obtained from the present frequency domain Volterra model describes the nonlinear transfer functions or Volterra Transfer Function (VTF) of the process in varied conditions of the patient and can be directly used in internal model control.
This work has been pursued as a part of thrust area research in smart instrumentation and control systems at the Electrical Engineering Department, BESU supported by the Special Assistance Programme of UGC, Govt. of India.
|||D. C. Howey. “The Treatment of Diabetes Mellitus”; Journal of Medical Pharmacology, pp. 1-9, Feb 2002.|
|||R. S. Parker, F. J. Doyle III, and N. A. Peppas. “A Model-based Algorithm for BG Control in Type 1 Diabetic Patients”; IEEE. Trans. on Biomedical Engg. Vol. 46, No.2, Feb 1999, pp. 148-157.|
|||R. S. Parker, F. J. Doyle III, J. H. Ward and N. A. Peppas. “Robust H∞ Glucose Control in Diabetes Using a Physiological Model”, Bioengineering, Food, and Natural Products, Vol. 46, No. 12, December 2000.|
|||A. Sutradhar and A. S. Chaudhuri; “Design and Analysis of an Optimally Convex Controller in Algebraic Framework for a Micro-insulin Dispenser System”, Journal of AMSE, France, Advances in Modeling & Analysis, Series C, Vol. 57 No. 1-2, 2002, pp. 1-14.|
|||A. Sutradhar and A. S. Chaudhuri; “Adaptive LQG/LTR Controller for Implantable Insulin Delivery System in Type-1 Diabetic Patient”; Proceedings of 3rd International Conference on System Identification and Control Problems, SICPRO’04, held at Institute of Control Sciences, Moscow, 28-30 January, 2004. pp. 1313-1328.|
|||R. Hovorka, V. Canonico, L. J. Chassin, U. Haueter, M. Massi- Benedetti, M. O. Federici, T., R. Pieber, H. C. Schaller, L. Schaupp, T. Vering, and M. E. Wilinska, "Nonlinear model predictive control of glucose concentration in subjects with type 1 diabetes," Physiological measurement, vol. 25, no. 4, Aug. 2004, pp. 905-920.|
|||I. Cochin and W. Cadwallender. Analysis and Design of Dynamic Systems, 3e, Addison-Wesley, ©1997.|
|||B. P. Kovatchev, M. Breton, C. D. Man, and C. Cobelli; “In Silico Preclinical Trials: A Proof of Concept in closed-loop control of diabetes”; Journal of Diabetes Science and Technology, Volume 3, Issue 1, January 2009, pp. 44-55.|
|||J. R. Guyton, R.O. Foster, J.S. Soeldner, M.H. Tan, C.B. Kahn, JL. Koncz, and R.E. Gleason, “A model of glucose-insulin homeostasis in man that incorporates the heterogeneous fast pool theory of pancreatic insulin release," Diabetes, vol. 27, 1978, pp. 1027-1042.|
|||J. T. Sorensen, “A Physiologic Model of Glucose Metabolism in Man and its Use to Design and Assess Improved Insulin Therapies for Diabetes”, Ph.D. thesis, Department of Chemical Engineering, MIT, 1985.|
|||E. D. Lehmann and T. Deutsch. “A Physiological Model of Glucose-insulin Interaction in Type-1 Diabetes Mellitus”; Jl. of Biomedical Engineering, Vol. 14, May 1992, pp. 235-242.|
|||C. D. Man, D. M. Raimondo, R. A. Rizza, and C. Cobelli; “GIM, Simulation Software of Meal Glucose–Insulin Model”; Journal of Diabetes Science and Technology, Volume 1, Issue 3, May 2007, pp. 323-330.|
|||A. Sutradhar and A. S. Chaudhuri; “Linear State-Space Model of Physiological Process in a Type-1 Diabetic Patient with Closed Loop Glucose Regulation”; Journal of AMSE, France, Advances in Modeling & Analysis, Series C, Vol. 66 No. 3, 2005, pp. 1-18.|
|||J. Li et. al. “Modeling the glucose–insulin regulatory system and ultradian insulin secretory oscillations with two explicit time delays”; Journal of Theoretical Biology 242, 2006 Elsevier; pp. 722-735.|
|||J. D. Chiu et. al. –“Direct Administration of Insulin Into Skeletal Muscle Reveals That the Transport of Insulin Across the Capillary Endothelium Limits the Time Course of Insulin to Activate Glucose Disposal” Diabetes, April 2008 vol. 57 no. 4, pp. 828-835.|
|||M. Sjostrand et. al. “Delayed Transcapillary Delivery of Insulin to Muscle Interstitial Fluid After Oral Glucose Load in Obese Subjects”: Diabetes, January 2005 vol. 54 no. 1, pp.152-157.|
|||B. Aussedat, M. Dupire-Angel, R. Gifford, J. C. Klein, G. S. Wilson, and G. Reach; “Interstitial glucose concentration and glycemia: implications for continuous subcutaneous glucose monitoring” American Journal of physiology Endocrinology and metabolism, April 1, 2000 vol. 278 no. 4 pp. E716-E728.|
|||R. S. Parker, F. J. Doyle III, and N. A. Peppas.”The intervenous route to blood glucose control”; IEEE Engineering in Medicine and Biology, Jan/Feb 2001; pp. 65-71.|
|||L. Ljung, System identification: Theory for the user, Prentice-Hall (Englewood Cliffs, NJ), 1987|
|||W. Zhao and H. Chen, Recursive identification for Hammerstein systems with ARX subsystem, IEEE Trans. Automat. Contr., 51(12), 2006, 1966-1974.|
|||A. Sutradhar and A. S. Chaudhuri; “On-line Identification of ARX model for Glucose-Insulin Interaction in a Diabetic Patient from Input-Output Data”; Journal of Institution of Engineers (India), Vol. ID-85, May 2004. pp. 1-7.|
|||A. Bhattacharjee, A. Sengupta and A. Sutradhar, “Nonparametric Modeling of Glucose-Insulin Process in IDDM Patient using Hammerstein-Wiener Model”; Proceedings of the 11th International Conference on Control, Automation, Robotics and Vision (ICARCV 2010), Singapore, 07-10 December 2010.|
|||D. T. Westwick and R. E. Kearney; “Identification of Nonlinear Physiological Systems”; -Wiley-Interscience, ©2003.|
|||G. Budura, C. Botoca, ‘Efficient Implementation and Performance Evaluation of the Second Order Volterra Filter Based on the MMD Approximation’, WSEAS Transactions on Circuits and Systems, Issue 3, Volume 7, March 2008; pp. 139-149.|
|||K. Sam Shanmugam, Mark T. Jong, “Identification of Nonlinear Systems in Frequency Domain” IEEE Transactions on Aerospace and Electronic Systems, Vol. AES-11, No.6 November 1975.|
|||G. Bicken, G. F. Carey and R.O. Stearman; “Frequency Domain Kernel Estimation for 2nd-order Volterra Models using Random Multi-tone Excitation”, VLSI Design; Taylor & Francis, Vol. 15(4), 2002, pp. 701-713.|
|||A. Bhattacharjee and A. Sutradhar, “Frequency Domain Hammerstein Model of Glucose-Insulin Process in IDDM Patient”; Proceedings of the International Conference on Systems in Medicine and Biology (ICSMB 2010), IIT Kharagpur, 2010.|
|||A. Bhattacharjee and A. Sutradhar, “Nonparametric Identification of Glucose-Insulin Process in IDDM Patient with Multi-meal Disturbance”; Journal of the Institution of Engineers (India): Series B, Vol. 93, No.4, pp. 237-246, Springer, Feb 2013.|
Dynamic Model of the type-1 diabetic patient
The differential mass balance equations governing the glucose-insulin process as per the compartmental model and meal model discussed by Parker , Guyton, Sorensen  and Lehmann  are given into two subgroups as follows.
A= auxiliary equation state (dimensionless).
B= fractional clearance (I, dimensionless; N, L/min).
G= glucose concentration (mg/dL).
I= insulin concentration (mU/L).
N= glucagon concentration (normalized, dimensionless).
Q= vascular plasma flow rate (L/min).
q= vascular blood flow rate (dL/min).
T= transcapillary diffusion time constant (min).
V = volume (L).
v = volume (dL).
= metabolic source or sink rate (mg/min or mU/min).
Glucose Model Sub- and Superscripts
A= hepatic artery
BU= brain uptake
C= capillary space
H= heart and lungs
HGP = hepatic glucose production
HGU = hepatic glucose uptake
IHGP= insulin effect on HGP
IHGU= insulin effect on HGU
IVI= intravenous insulin infusion
K = kidney
KC= kidney clearance
KE= kidney excretion
LC= liver clearance
NHGP = glucagon effect on HGP
P= periphery (muscle/adipose tissue).
PC= peripheral clearance
PGU= peripheral glucose uptake
PIR = pancreatic insulin release
PNC= pancreatic glucagon clearance
PNR= pancreatic glucagon release (normalized).
RBCU= red blood cell uptake
S= gut (stomach/intestine).
SIA = insulin absorption into blood stream from subcutaneous depot
SU= gut uptake
T= tissue space
Glucose sub-model mass balance equations:
The metabolic source and sink terms () in mg/min and the uncertainty equation given in section V with A, B, C, D having usual significances are taken from . represents the rate of glucose into the gut compartment of diabetic patient where is the input pulse signal.
Insulin sub-model mass balance equations:
The metabolic source and sink terms () in mU/min. overall type-1 diabetic patient has been characterized by 19 differential equations, with eleven describing glucose dynamics, seven for insulin dynamics, and a single compartment for glucagon . Model parameters are as shown in Table A1.