Glucose optimal control system in diabetes treatment

Applied Mathematics and Computation 209 (2009) 19–30

Contents lists available at ScienceDirect

Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc

Glucose optimal control system in diabetes treatment Irma Y. Sánchez Chávez a,*, Rubén Morales-Menéndez b, Sergio O. Martínez Chapa c a b c

Mechatronics and Automation Department, Tecnológico de Monterrey, Eugenio Garza Sada 2501 Sur, 64849 Monterrey NL, Mexico Center for Innovation in Design and Technology, Tecnológico de Monterrey, Eugenio Garza Sada 2501 Sur, 64849 Monterrey NL, Mexico Electrical Engineering Department, Tecnológico de Monterrey, Eugenio Garza Sada 2501 Sur, 64849 Monterrey NL, Mexico

a r t i c l e

i n f o

Keywords: Optimization Glucose optimal control Closed loop drug delivery

a b s t r a c t Optimal control allows the incorporation of functional constraints and requirements as a departure point for the design process. A control system for optimal insulin delivery in a type I diabetic patient is presented based on the linear quadratic control problem theory. The glucose–insulin dynamics is first represented by a linear model whose state variables are the glucose and the insulin concentrations in the blood. These variables allow the formulation of an appropriate cost function for a diabetes treatment in terms of the deviation from the normal glucose level and the dosage of exogenous insulin. The optimal control law is computed from this cost function under the servocontrol and regulatory approaches. A Monte Carlo simulation shows the superior robustness of the regulatory control design. Further evaluation of the regulatory controller is realized with a high order non-linear human glucose–insulin model. The control system performance can be improved by adjusting the weighting factors of the optimization problem according to the patient’s needs. The velocity of the correction of a glucose level deviation can be increased with a higher value for its weighting factor, while an augment of the weighting factor for insulin supply may be necessary to prevent saturation of the controller output and oscillation of the blood glucose concentration. Results exemplify the suitability of the optimization and regulatory approaches to produce a practical control algorithm for biomedical engineering problems. Ó 2008 Elsevier Inc. All rights reserved.

1. Introduction Computational tools have enabled the simulation of complex processes for analysis and control. Automatic control techniques can improve the functionality of biomedical systems, including disease treatments. Diabetes mellitus has been one of the main causes of death in Mexico in the last years, and this fact inserts into a worldwide reality of an increasing number of cases. Juvenile-onset diabetes is characterized by the defective or null capacity to produce insulin, necessary for the absorption of glucose in the cells. This problem itself determines a medical treatment based on the dosing of exogenous insulin, which is done frequently in a fixed prescribed amount without continuous monitoring or feedback information. Feedback control systems have the potential for a precise blood glucose regulation in a diabetic patient. Insulin infusion can be done subcutaneously and intravenously. The subcutaneous route is easier and safer to manage, which is an advantage for closed systems implementation since it is the route used in traditional open loop diabetes treatment. The intravenous route avoids time delays to reach blood stream and to produce body response, which is convenient for continuous closed loop performance. Both insulin delivery types have been considered in closed loop treatment systems [1,2]. Lispro insulin may combine the advantages of both routes because of its fast absorption after subcutaneous injection. * Corresponding author. E-mail address: [email protected] (I.Y. Sánchez Chávez). 0096-3003/$ - see front matter Ó 2008 Elsevier Inc. All rights reserved. doi:10.1016/j.amc.2008.06.030

20

I.Y. Sánchez Chávez et al. / Applied Mathematics and Computation 209 (2009) 19–30

An open loop treatment system may be seen as a partially closed loop system because medical prescription of insulin is based on home glucose monitoring among other information of the patient [1]. Decision support systems have been designed for diabetes management for this kind of therapy systems. Controllers defined with proportional, integral and derivative terms of the difference between the desired and the actual process variable values (PID controllers) have been used for blood glucose control. A PID based on a sliding scale approach, tested in patients in intensive care unit, has been reported [3]. A proportional-derivative or PD controller has been derived with a pole assignment strategy and tested in patients [1]. Robustness of closed loop systems has been achieved by adaptive and predictive mechanisms to account for sparse glucose measurement [4] and time variations of the glucose–insulin process. The model predictive control algorithm implements a self-tuning controller that has been studied by simulation [5] and tested in vivo [6]. In this work, a Linear Quadratic Control (LQC) problem is stated for the design of an optimal controller. A quadratic performance criterion is expressed in terms of the blood glucose deviation from its normal value and the amount of exogenous insulin to supply. The controller design, therefore, aims to improve the condition of the patient with the minimum dose of medicament. The formulation approach of the performance criterion for an optimal system may lead to tracking problems [2]. Servocontrol and regulatory approaches are discussed for automatic insulin delivery. The glucose–insulin physiologic process to control is first represented with a linear second order model [7], then a nonlinear 19th order model [8] is used for a more detailed description and closer proximity to the reality and complexity of this process. The system is closed by coupling the optimal controller with the high order, non-linear physiological process. The effect of the weighting factors of the optimization problem in the closed loop response is also analyzed. The paper is organized as follows: Section 2 presents the mathematical models for the glucose–insulin dynamics. In Section 3, a blood glucose level optimal controller is designed and a state observer is suggested. The simulation results of the proposed closed loop treatment system are presented in Section 4. Finally, in Section 5, the paper is concluded. 2. Mathematical models of glucose metabolism 2.1. Ackerman’s glucose–insulin model Research on glucose control in diabetic patients depends on the development of mathematical models for the glucose– insulin dynamics in the human body. Compartmental models are based on mass balances between different organs or compartments. The Ackerman’s model is widely known because of its simplicity, since it considers one compartment that represents the global glucose–insulin dynamics in the human body [7]. This model is based on the glucose tolerance test where the basal level is disturbed by the intake of glucose. The non-linear interaction between glucose and insulin is described by

G0 ðtÞ ¼ f1 ðG; HÞ þ pðtÞ;

ð1aÞ

H0 ðtÞ ¼ f2 ðG; HÞ þ uðtÞ;

ð1bÞ

with G(t = 0) = G0, H(t = 0) = H0, p(t = 0) = 0 and u(t = 0) = 0, where G(t) is the glucose level, H(t) is the hormone level, p(t) is the external glucose supply rate and u(t) is the insulin infusion rate at time t. The linearization of the model is necessary in order to obtain a closed solution to the optimal control problem. Considering the deviation variables g(t) = G(t)  G0 and h(t) = H(t)  H0, and applying a standard linearization procedure, the linear representation of the model is

g 0 ðtÞ ¼ m1 gðtÞ  m2 hðtÞ þ pðtÞ; 0

h ðtÞ ¼ m4 gðtÞ  m3 hðtÞ þ uðtÞ:

ð2aÞ ð2bÞ

The terms containing m1 and m3 account for self-removal of glucose and insulin, respectively, therefore m1 and m3 are positive. The terms with m2 and m4 represent the reduction of the glucose level with insulin and the increment of blood insulin with glucose, respectively, which implies that m2 is positive and m4 is greater than or equal to zero. The parameters m1, m2, m3 and m4 have been obtained from experimental data [9]. In the case of a type I diabetic patient, m4 = 0. The reported values in min1 for the rest of the parameters are m1 = 0.0009, m2 = 0.0031 and m3 = 0.0415. The term p(t) is not considered for the controller design. The final state-space model is given by

x0 ðtÞ ¼ AxðtÞ þ BuðtÞ;

xð0Þ ¼ 0;

yðtÞ ¼ CxðtÞ þ DuðtÞ;

ð3aÞ ð3bÞ

where x(t) = [x1(t) x2(t)]T = [g(t) h(t)]T, A = [m1 m2; 0 m3], B = [0 1]T, C = [1 0] and D = 0. 2.2. Sorensen’s glucose–insulin model The compartmental technique applied to the main organs of the human body involved in the glucose–insulin dynamics, the consideration of convection and diffusion transport mechanisms and the representation of the underlying kinetics have

I.Y. Sánchez Chávez et al. / Applied Mathematics and Computation 209 (2009) 19–30

21

led to the non-linear 19th order model contributed by Sorensen [8]. This model has been the base for research on meal and exercise effects on glucose–insulin metabolism and investigation on blood glucose control. The present work takes the Sorensen’s model with the inclusion of meal disturbance modelling and the corresponding physiologic parameters available in the literature [10]. Fig. 1 shows a schematic representation of the Sorensen’s model. The application of the solution of the LQC problem requires a linear approximation of the Sorensen model. In order to fit the structure of the linear Ackerman model (3), the transition matrix A is obtained by the evaluation of the Jacobian matrix:

 A¼

m1

m2

m4

m3



 ¼

of1 =ox1

of1 =ox2

of2 =ox1

of2 =ox2

 ð4Þ

with m4 = of2/ox1 = 0, where f1 and f2 are the first time derivatives of glucose (x1) and insulin (x2) concentrations, respectively, in the heart and lungs compartment. This linearization procedure preserves meaningful parameters and states as well as mathematical precision, since the simplified model considers explicitly the main physiologic variables in diabetes. 3. Optimal controller 3.1. Linear quadratic control problem The LQC problem consists of determining a control law u(t) to minimize the cost function given by the next equation:

JðuÞ ¼

1 T 1 e ðt f ÞSeðt f Þ þ 2 2

Z

tf

½eT ðtÞQðtÞeðtÞ þ uT ðtÞRðtÞuðtÞdt;

ð5Þ

t0

where u(t) appears quadratically; e(t) is the error vector, S is a constant matrix; Q and R may vary with time; Q, R and S are symmetric matrices; S and Q are positive semidefinite, and R is positive definite. The control law u(t) is the input of the state-space model:

x0 ðtÞ ¼ AðtÞxðtÞ þ BðtÞuðtÞ;

xðt 0 Þ ¼ x0 ;

t 2 bt0 ; t f c;

ð6Þ

where u(t) is included linearly, and the state transition matrix A(t) and the input matrix B(t) may be time variant in general. An optimal control law u*(t) is assumed to exist for this problem in [t0, tf]. An optimal trajectory x*(t) is associated with u*(t). The control signal and the state vector can be expressed as

uðtÞ ¼ u ðtÞ þ etðtÞ;

ð7aÞ

xðtÞ ¼ x ðtÞ þ eyðtÞ;

ð7bÞ

where e is a small positive number and t(t) is arbitrary. The optimal control law is obtained when e = 0, which makes dJ(e)/de = 0 [11,12]. The solution consists of the control law:

uðtÞ ¼ RðtÞ1 BðtÞT kðtÞ;

ð8Þ

where k(t) is the coestate vector defined as

Brain BGC mesurement Venous blood Heart/Lungs Insulin delivery

Hepatic Artery

Liver

Gut Portal Vein

Meal disturbance

Kidney Periphery Fig. 1. The Sorensen’s compartmental model for the glucose–insulin metabolism (BGC: blood glucose concentration).

22

I.Y. Sánchez Chávez et al. / Applied Mathematics and Computation 209 (2009) 19–30

kðtÞ ¼ PðtÞxðtÞ þ lðtÞ;

ð9Þ

where P is a symmetric positive semidefinite matrix and l is a column vector and the dynamic behaviour of the coestate vector is given by

k0 ðtÞ ¼ Q ðtÞeðtÞ  AT ðtÞkðtÞ:

ð10Þ

The existence of an independent term in the control law depends on the exploited design approach in the formulation of the linear quadratic problem [13]. 3.2. Servocontrol design Under the servocontrol approach, let xd be the desired glucose change in a diabetic person (xd = Gd  G0). Deviation of blood glucose level from its desired value (x1  xd = G  Gd) and insulin infusion must be minimized. The objective function is

JðuÞ ¼

Z

1

½gðx1 ðtÞ  xd Þ2 þ qu2 ðtÞdt;

ð11Þ

0

where g and q are positive weighting factors. By comparing equations (5) and (11), it can be identified that S ¼ ½ 0 0; 0 0 , Q ¼ ½ 2g 0; 0 0 , R = 2q, t0 = 0 and tf = 1. The open loop system is represented by a state space model with constant parameters. The closed loop dynamics achieved with the optimal control law is obtained by combining equations (6) and (8):

x0 ðtÞ ¼ AxðtÞ  BR1 BT kðtÞ:

ð12Þ

From Eqs. (9) and (10), the first derivative of the coestate vector is

Px0 ¼ Q ðtÞ½xðtÞ  xD   AT PxðtÞ  AT l;

ð13Þ

where xD(1) = xd. Replacing x0 (t) with Eq. (12) in (13) gives

ðPA þ AT P  PBR1 BT P þ Q ÞxðtÞ ¼ ðPBR1 BT  AT Þl þ QxD :

ð14Þ

The matrix P is chosen to satisfy the Riccati equation:

PA þ AT P  PBR1 BT P þ Q ¼ 0:

ð15Þ

Then, from Eq. (14), l can be calculated (l – 0). The optimal control law is

uðtÞ ¼ K c x þ K; where K c ¼ b K c1

K c2 ¼

ð16Þ

K c2 c and the parameters are calculated as

rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi m21 þ m23 þ 2 m21 m23 þ m22 q1 g  ðm1 þ m3 Þ;

ð17Þ

K c1 ¼ ðK 2c2 þ 2m3 K c2 Þ=2m2 ;

ð18Þ

K ¼ m2 xd q1 ðm21 m23 þ m22 q1 gÞ0:5 :

ð19Þ

With g = 1, q = 10 min2 and xd = 200 mg/dl (according to the simulation problem presented in Section 4) and the Ackerman model parameters of Section 2.1, the resulting parameters of the control law are Kc1 = 0.2991 min1, Kc2 = 0.0183 min1 and K = 63.1997 mg/dl/min. 3.3. Regulatory design For the regulatory approach, a normal glucose level is referred as the initial steady state. Any deviation from this value is a disturbance. Therefore, x1 is the deviation of glucose level from the desired value (Gd = G0; xd = 0; x1 = G  Gd) and the performance criteria is expressed as

JðuÞ ¼

Z 0

1

½gx21 ðtÞ þ qu2 ðtÞdt:

ð20Þ

Considering the linear state space process model, the performance criteria represented by Eq. (20) and the application of Eqs. (12)–(15), l = 0 and the control law is u(t) = Kcx(t). With g = 1, q = 10 min2 and the Ackerman model parameters of Section 2.1, Kc = [0.2991 0.0183].

I.Y. Sánchez Chávez et al. / Applied Mathematics and Computation 209 (2009) 19–30

23

3.4. State estimation The control law in both approaches is a function of the two states of the system, glucose (x1) and insulin (x2) concentrations. The glucose concentration measurement is supposed to be available while the insulin concentration needs to be estimated. The design of the control law and the design of the state observer are independent. The observer model is given by

~x0 ðtÞ ¼ A~xðtÞ þ BuðtÞ þ K e ½yðtÞ  C ~xðtÞ;

ð21Þ

where yðtÞ  C ~ xðtÞ is the observation error. The Ackerman’s model is an observable system. The state observer is designed according to the classical control theory using the parameters given in Section 2.1. The gain vector Ke = [0.2021 2.1023]T produces a faster observer response (natural frequency xn = 0.122 rad/min and damping ratio f = 1), with respect to the closed loop behaviour (xn = 0.031 rad/min and damping ratio f = 0.97). If the physiologic process is represented by the Sorensen non-linear model, the best linear approximation may vary with time, which causes ‘‘process noise” [14]. Additionally, noisy sensor measurements can be considered to justify the use of a Kalman filter for the sate estimation. 4. Computer simulation The simulation test consists of a situation of hyperglycaemia in a type I diabetic patient with an initial glucose level of 300 mg/dl. The closed loop system should reach a desired steady state level of 100 mg/dl. The sensor and actuator are assumed to be ideal systems and the open loop dynamics is supposed to be dominated by the model of the patient or physiologic process to be controlled. 4.1. Closed loop with a linear process model The Ackerman linearized model for the glucose–insulin process can be used for the discussion of control issues, which can be solved before using a broader model. The purpose of the simulation of the linear closed loop system is to evaluate the controller design approach in the correction of hyperglycaemia. The recovery of a normal glucose level of 100 mg/dl with the elimination of the cause of the hyperglycaemic condition and without supply of exogenous insulin is shown in Fig. 2. The control loop is closed with the classical state observer and the optimal controller as illustrated in Fig. 3. This configuration allows the comparison with reported results for the stated simulation problem [15]. Under a servocontrol approach, the initial deviation value is zero assuming an initial stable state in the hyperglycaemic condition (300 mg/dl). The glucose change to be achieved is 200 mg/dl. Simulation results are shown in Fig. 4. The transitory period elapses 4 h, a steady state error of 0.3 mg/dl is detected and the cost function grows indefinitely. For the regulatory approach, the reference steady state is the normal condition of glucose level at 100 mg/dl, so the initial deviation of 200 mg/dl is considered a disturbance to the closed loop system. No offset error is obtained and the cost function converges, as Fig. 5 shows. A linearized and low order glucose–insulin model may cause uncertainty as non-linear high order models do. Although a more complete model may be suitable to manage characterized uncertainty [10], all effects may not be represented, which justifies a random variation of model parameters to analyze controller robustness.

Fig. 2. Open loop response using the Ackerman linearized model for the glucose–insulin metabolism.

24

I.Y. Sánchez Chávez et al. / Applied Mathematics and Computation 209 (2009) 19–30

Optimal controller (LQC)

BIC estimation

BGC Glucose-insulin measurement metabolism in a diabetic patient Ackerman linearized model

Insulin delivery

State observer

Fig. 3. Closed loop system for the control of the glucose–insulin metabolism given by the Ackerman linearized model (BGC: blood glucose concentration; BIC: blood insulin concentration).

Fig. 4. Servocontrol results with a linearized model for glucose–insulin metabolism: (a) transient response; (b) final steady state; (c) cost function.

A Monte Carlo simulation proves the robustness of the control system. The parameters of the Ackerman model are changed following a stochastic process, while the controller parameters are not altered. To illustrate the variability in the performance of the closed loop system before changes in process dynamics, Fig. 6 presents box and whisker plots for the cost function and the offset of the closed loop response under different percentages of variation. These graphs were computed by 30 independent runs of the simulation test. Beyond 60% parameter variation, the servocontrol system shows unacceptable performance. The regulatory design performs with no significant difference with more than 50% variation.

I.Y. Sánchez Chávez et al. / Applied Mathematics and Computation 209 (2009) 19–30

25

Fig. 5. Regulatory results with a linearized model for glucose–insulin metabolism: (a) transient response; (b) final steady state; (c) cost function.

4.2. Closed loop with a non-linear process model A second closed loop system is proposed using the Sorensen’s model to simulate the glucose–insulin metabolism in a diabetic patient. This model resembles the non-linear nature of the real process. The model response converges to a normal steady state glucose concentration of 119.5 mg/dl with no infusion of insulin. Therefore the simulation problem is slightly modified by the consideration of the desired final glucose concentration of 119.5 mg/dl. The hyperglycaemic condition is produced by the simulation of a glucose meal input to the gut compartment of the Sorensen’s model. The glucose meal input is increased exponentially from zero and kept at 408.7 mg/min to stabilize the arterial glucose concentration in 299.65 mg/ dl. Then, the glucose meal input is suspended and the metabolism of the diabetic patient takes about 8 h to reach a steady glucose concentration of 119.5 mg/dl without insulin delivery, as shown in Fig. 7. The control system includes a linearization module, a Kalman filter and the optimal controller, as indicated in Fig. 8. The linearization module is used for the calculation of the parameters of the control law and the Kalman estimator gain vector. The Kalman filter is needed to estimate the blood insulin concentration. The control approach to apply is the regulatory one, since it produces higher performance with a simpler configuration (Section 4.1). The simulation diagram of Fig. 8 helps to evaluate the adequacy of the control law structure and regulatory approach for a non-linear process in the presence of noise. The weighting factors g and q of the cost function are adjusted to avoid oscillation. The normal glucose level can be restored after 250 min, according to the simulation results shown in Fig. 9. The closed loop system reduces the open loop settling time (Fig. 7), but the difference is less drastic than in the case of the system of Fig. 3, as transient responses in Figs. 2 and 5 show, since the controller is not based on the same model used for the glucose–insulin process simulation.

26

I.Y. Sánchez Chávez et al. / Applied Mathematics and Computation 209 (2009) 19–30

Fig. 6. Robustness tests: (a) servocontrol approach; (b) regulatory approach. Cost function is J given by Eqs. (11) or (20) in (mg/dl)2 min  106; offset is reported in mg/dl. The boxes have lines at the lower quartile, median, and upper quartile values. The whiskers are lines extending from each end of a box to show the extent of the rest of the data. The boxes are notched. Notches represent a robust estimate of the uncertainty about the medians for box-to-box comparison. Outliers are data values beyond the ends of the whiskers; outliers are shown as circles.

Fig. 7. Open loop response using the Sorensen non-linear model for glucose–insulin dynamics.

Optimal controller (LQC)

Insulin delivery

Glucose-insulin metabolism in a diabetic patient

BGC measurement

Sorensen nonlinear model BIC estimation

Kalman filter

Ackerman linear model parameters

Linear model fitting

Fig. 8. Closed loop system for the control of the glucose–insulin metabolism represented with the Sorensen non-linear model.

I.Y. Sánchez Chávez et al. / Applied Mathematics and Computation 209 (2009) 19–30

27

Fig. 9. Regulatory approach results with the Sorensen non-linear model for the glucose–insulin dynamics: (a) transient response; (b) final steady state; (c) cost function.

The parameters g and q of the cost function (20) are varied to show how they affect the system dynamics in Fig. 10. The magnitude of the initial slope and the undershoot of the closed loop response increase with higher g and lower q. Fig. 10 also shows the corresponding effect of the weighting factors on insulin delivery (bottom plots). The order of magnitude of q needs to be very high in order to prevent saturation and an ON–OFF controller behavior. The values for the weighting factor q are very different from the value used when the linear process model represents the glucose–insulin metabolism due to the distinct dynamics and different units required for the metabolic source of insulin in the Sorensen model. 4.3. Summary of results The simulation of a closed loop system for diabetes treatment considers a glucose–insulin process model and a controller. The low order and linearized Ackerman model provides a practical basis for the evaluation of the servocontrol and regulatory controller design approaches. Simulation results with the Ackerman model and different design approaches show the superiority of the regulatory controller, since it avoids a steady state error in the glucose level (Figs. 4 and 5) and realizes corrections and maintains stability in spite of dynamic variations of the physiological process, which are examined by means of a Monte Carlo simulation (Fig. 6). The advantageous regulatory design is applied to control the glucose–insulin metabolism represented by the Sorensen high order, non-linear model. The closed loop system is refined by considering a more detailed description of the process, a Kalman filter to account for the noise effect in the estimation of the blood insulin concentration, and an adaptive mecha-

28

I.Y. Sánchez Chávez et al. / Applied Mathematics and Computation 209 (2009) 19–30

Fig. 10. Effect of weighting factors on closed loop response. (a) The values for g are 1, 1.5, 2 and 2.5, with q = 3.3  1011 min2. (b) The factor q is varied as 1, 2, 3 and 4  1011 min2 with g = 1.

nism to update the controller parameters (Fig. 8). The added complexity responds to basic requirements for a real application of the proposed system. The results are summarized in Table 1. The possibility for external adjustment by the patient or the physician, as a safety feature or extraordinary measure, is also considered by means of the weighting factors of the performance function. A low weighting factor for the glucose level deviation and a high weighting factor for the exogenous insulin lead to a slower correction of the blood glucose concentra-

I.Y. Sánchez Chávez et al. / Applied Mathematics and Computation 209 (2009) 19–30

29

Table 1 Simulation results for the correction of a hyperglycaemic condition of 300 mg/dl Physiological model Open loop settling time Control approach Closed loop settling time Closed loop offset

Ackerman 3500 min Servocontrol 230 min Variable

Ackerman 3500 min Regulatory 200 min Null

Sorensen 450 min Regulatory 250 min Null

tion with more restricted supply of insulin (Fig. 10). Appropriate weighting factors would give a faster response without an oscillation that could imply a hypoglycaemia risk. 5. Conclusions The stated simulation test with the system configuration of Fig. 3 has been presented in the literature [15]. Glucose levels at 1 and 2 h with the servocontrol approach are very similar to reported values. Glucose levels of 185.98 and 119.86 mg/dl are obtained while previous works report 183 and 119 mg/dl at these times, respectively. For the regulatory approach, glucose concentrations of 166.8 and 112.57 mg/dl are reached at 60 and 120 min. These results differ more from the reported results, however, the performance is more satisfactory. The structure of the control law obtained with the LQC technique corresponds to a proportional controller, which causes an offset except in the steady state condition of the patient used as reference to define deviation values. If the reference steady state of the patient is the state of normal blood glucose level, the control action proportional to the states can reach the desired glucose level from any abnormal condition. Monte Carlo simulations prove the higher robustness for the regulatory control scheme. The closed loop implemented with the process dynamics given by the Sorensen’s model under a regulatory approach (with g = 1 and q = 3.3  1011 min2) achieves 60.87% and 93.06% of the total desired change at 60 and 120 min, with no steady state deviation. Reported results show a change in blood glucose concentration of 58.5% and 90.5% at the same times. Thus, the comparison shows agreement with prior researches and validates the present approach for optimal glucose control. Although a closed loop treatment system is designed to perform without human intervention, the user (patient or physician) could make adjustments by modifying the weighting factors g for the glucose concentration deviation and q for the insulin delivery or controller action. Since the inverse of the weighting factor for the insulin delivery and the weighting factor for the glucose level deviation have similar effects, and the product of both appears in the calculations of the controller gain vector Kc, a single means can be proposed for tuning the system, either q or g or the product q1g itself. A unique parameter for adjustment would ease the interaction of the user with the system. The results based on the interaction of the optimal controller with a highly non-linear plant show the robustness of the controller and the adequacy of the design approach for a particular biomedical problem. However, the proposed system lacks of a feedforward or anticipation mechanism that is present even in the traditional open loop treatment based on insulin injections previous to the ingestion of food, for example. The presented closed loop system intends to emulate the response of a healthy body to the blood glucose level resultant from the activity of the patient, trying not to alter his daily routine. Practical limitations arise regarding the implementation of the closed loop system. The system relies on a continuous meter of glucose levels and an analogous insulin pump. Although both types of devices have been studied and developed, their continuous operation still has to be improved for a smooth functioning [16]. The models of a glucose sensor and an actuator to deliver insulin should be considered specifically for an evaluation of the closed loop treatment from the perspective of the comfort of the patient, as well as from the dynamics point of view in future work. In addition to the instrumentation challenges, unavailable measurements and the complexity of physiological processes demand reliable state estimation which needs noise management and on line identification. Further investigation could be oriented to the use of unscented particle filters for non-linear processes to improve blood insulin state estimation [17]. The tests for the closed control system should include the regulation during prolonged disturbances such as exercising or feeding which may require the extension of the physiologic model for simulation analysis. Continuous feedback control systems for disease treatment are potentially bloodless, painless and more precise therapies. Their successful application in biomedicine requires the development of microsystems with reliable sensors and actuators and embedded control algorithms. Modelling and computer simulation assist the integration of sensing, actuating and controlling components and the evaluation of their interaction with the process or medical situation. Therefore, interdisciplinary collaboration is needed for research and development of biomedical systems. Acknowledgements The authors thank Dr. Graciano Dieck Assad, Dr. Frantz Bouchereau Lara and Prof. Óscar Miranda Domínguez for their valuable comments and the Consejo de Ciencia y Tecnología del Estado de Nuevo León for financial support.

30

I.Y. Sánchez Chávez et al. / Applied Mathematics and Computation 209 (2009) 19–30

References [1] R. Bellazi, G. Nucci, C. Cobelli, The subcutaneous route to insulin-dependent diabetes therapy, IEEE Eng. Med. Biol. 20 (1) (2001) 54–64. [2] R. Parker, F. Doyle III, N. Peppas, The intravenous route to blood glucose control, IEEE Eng. Med. Biol. (2001) 65–73. [3] F. Chee, T. Fernando, P. Van Heeden, Expert PID control system for blood glucose control in critically ill patients, IEEE Trans. Inform. Tech. Biomed. 7 (4) (2003) 419–425. [4] E. Woodruff, S. Gulaya, R. Northrop, The closed-loop regulation of blood glucose in diabetes, in: Proceedings of the 14th Annual Northeast Bioengineering Conference, 1988, pp. 54–57. [5] R. Parker, F. Doyle III, N. Peppas, Model-based algorithm for blood glucose control in type i diabetic patients, IEEE Trans. Biomed. Eng. 46 (2) (1999) 148–157. [6] R. Dudde, T. Vering, Advanced insulin infusion using a control loop (ADICOL) concept and realization of a control-loop application for the automated delivery of insulin, in: Proceedings of the Fourth International IEEE EMBS Special Topic Conf. on Info. Tech. Appl. in Biomed., 2003, pp. 280–282. [7] E. Ackerman, L. Gatewood, J. Rosevear, G. Molnar, Model studies of blood-glucose regulation, Bull. Mathem. Biophys. 27 (Suppl.) (1965). [8] J. 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. [9] T. Yipintsol, L. Gatewood, E. Ackerman, P. Spivak, G. Molnar, J. Rosevear, F. Service, Mathematical analysis of blood glucose and plasma insulin responses to insulin infusion in healthy and diabetic subjects, Comput. Biol. Med. 3 (1973) 71–78. [10] R. Parker, F. Doyle III, J. Ward, N. Peppas, Robust H1 glucose control in diabetes using a physiological model, AIChE J. 46 (12) (2000) 2537–2549. [11] D. Naidu, Optimal Control Systems, CRC Press, Boca Raton, FL, 2003. [12] R. Vinter, Optimal Control, Birkhöuser, Boston, MA, 2000. [13] I.Y. Sánchez Chávez, R. Morales-Menéndez, S.O. Martínez Chapa, Linear quadratic control problem in biomedical engineering, in: L. Puigjaner, A. Espuña (Eds.), Computer Aided Chemical Engineering, vol. 20: European Symposium on Computer Aided Process Engineering-15, Elsevier, 2005, pp. 1195–1200. [14] S. Haykin, Adaptive Filter Theory, Prentice Hall, Upper Saddle River, NJ, 2002. [15] M. Kikuchi, N. Machiyama, N. Kabei, A. Yamada, Y. Sakurai, Homeostat to control blood glucose level, in: Int. Symp. Med. Inf. Syst., 1978, pp. 541–545. [16] Diabetes Technology Meeting, Sixth Annual Meeting, Chair D. Klonoff, Atlanta, GA, 2006. [17] R. Morales-Menéndez, J. Mutch, N. de Freitas, D. Poole, F. Guedea-Elizalde, Dynamic modelling and control of industrial processes with particle filtering algorithms, in: Barbosa-Póvoa, H. Matos (Eds.) ESCAPE-14, Lisbon, Portugal, 2004, pp. 721–726.