Systematically In Silico Comparison of Unihormonal and Bihormonal Artificial Pancreas Systems

Automated closed-loop control of blood glucose concentration is a daily challenge for type 1 diabetes mellitus, where insulin and glucagon are two critical hormones for glucose regulation. According to whether glucagon is included, all artificial pancreas (AP) systems can be divided into two types: unihormonal AP (infuse only insulin) and bihormonal AP (infuse both insulin and glucagon). Even though the bihormonal AP is widely considered a promising direction, related studies are very scarce due to this system's short research history. More importantly, there are few studies to compare these two kinds of AP systems fairly and systematically. In this paper, two switching rules, P-type and PD-type, were proposed to design the logic of orchestrates switching between insulin and glucagon subsystems, where the delivery rates of both insulin and glucagon were designed by using IMC-PID method. These proposed algorithms have been compared with an optimal unihormonal system on virtual type 1 diabetic subjects. The in silico results demonstrate that the proposed bihormonal AP systems have outstanding superiorities in reducing the risk of hypoglycemia, smoothing the glucose level, and robustness with respect to insulin/glucagon sensitivity variations, compared with the optimal unihormonal AP system.


Introduction
According to a prediction produced by the International Diabetes Federation (IDF) in 2009, there will be approximately 285 million diabetic subjects by 2010 and about 439 million diabetic subjects by 2030 [1,2]. Diabetes has become a big issue that seriously harms human health. To regulate the blood glucose level [3,4], exogenous insulin infusion is one of our main choices. In the recent decades, there has been a great deal of interest in developing a closed loop system that embodies an automated insulin delivery system [5] for regulating blood glucose (BG) in type 1 diabetes. Unlike intravenous (IV) administration, however, the subcutaneous (SC) insulin infusion presents a control challenge owing to a finite absorption rate of insulin analogs. This challenge may induce excessive drug accumulation in the subcutaneous tissue and then even hypoglycemia (blood glucose lower than 70 mg/dL).
Short-term hypoglycemia may cause more damage to patients' brain function and even cause death [6]. Therefore, minimizing or avoiding hypoglycemia is a daily challenge for diabetic patients. There are mainly three ways to rescue for hypoglycemia: mathematically minimizing the SC accumulation, suspending the insulin delivery (termed as prediction/suspending solution) [7,8], and using glucagon as a counterregulatory agent in closed-loop system. The first method to handle hypoglycemia is to mathematically minimize the SC accumulation, for example, MPC [9][10][11]. We can design an objective function to optimally regulate the BG close to a target within the normoglycemic range and in the meantime to minimize the SC accumulation of insulin. Another strategy to avert hypoglycemia is a combination of glucose prediction [12] and insulin infusion suspending. The third solution is adding glucagon infusion [13].
Recently, El-Khatib et al. proposed a new scheme for glucose management [14]: dual subcutaneous infusion of insulin and glucagon, where a model prediction control (MPC) algorithm [15] was used to design the SC administration of insulin and a proportional-derivative control [16] was used to govern the glucagon infusion. Ward and his coworkers used the glucagon to prevent hypoglycemia in type 1 diabetes [17], where the fading memory proportional  Figure 1: The block diagram of simulation models. The simulation models for a virtual subject consist of a glucose subsystem, an insulin subsystem, a meal subsystem, and a glucagon subsystem. derivative (FMPD) algorithm was used to design the subcutaneous insulin and glucagon infusion rates. Several clinical trials [18,19] have been finished and they demonstrate the feasibility and safety of BG control by a bihormonal artificial pancreas. However, in the literature [14,[17][18][19][20], the delivery rates for insulin and/or glucagon were determined separately. Obviously, separate design cannot exploit the full benefit of using bihormone and may result in hormone waste and even blood glucose fluctuation.
In this paper, a new method based on switching control theory [21,22] and IMC-PID [23] controller was proposed to orchestrate switching between insulin and glucagon subsystems. Two kinds of switching rules were designed and compared in this paper [24]: P-type switching rule and PDtype switching rule. Compared with the reported bihormonal artificial pancreas systems [14,[17][18][19][20], the proposed algorithms have the following advantages. First, the proposed algorithms can coordinate the insulin and glucagon delivery simultaneously, such that more hypoglycemia can be minimized or avoided and more glucose concentrations can be kept within the normoglycemic range. Second, unreasonable operation, for example, delivering two hormones at the same time, was avoided in the new systems; hence, hormone waste and BG fluctuation can be minimized. Third, the transient phase (period during which neither insulin nor glucagon is delivered) can be further maximized, and the dosage of two hormones can be minimized.
Even though there are a number of reported works on clinical testing of bihormonal AP system [14,[17][18][19][20], its superiority compared with the uni-hormonal AP system is not fully proved and widely accepted. For example, in a comment on the bihormonal therapy [25], it is stated that "in terms of glucose control, the results were quite similar to uni-hormonal systems. " It is further suggested in [25] that "the efficacy of glucagon should be tested compared to an insulin-only system". Obviously, fair, systematic, and strict comparisons of the bihormonal and uni-hormonal AP systems are very important and valuable. However, strictly fair comparison of two therapies is impossible to carry out in clinic due to inevitable disturbances. Because any disturbance can be avoided in computer simulation, in silico test provides a possibility to fairly, systematically, and strictly compare of the bihormonal and uni-hormonal AP systems. In the in silico trials, various therapies can be tested in the same scenario.
Because IMC-PID is used to design the insulin and glucagon delivery rates in the proposed bihormonal AP systems, IMC-PID is also implemented in the uni-hormonal benchmark AP system. As introduced in the previous paragraph, adding prediction/suspending term can enhance the performance of the uni-hormonal AP system. There have been a number of reported studies for hypoglycemia prediction [7,8]. To achieve the performance limitation of the unihormonal AP system with prediction/suspending term, the perfect prediction/suspending solution is included in the unihormonal benchmark AP system; hence, this system can be considered as an optimal uni-hormonal AP system. To the authors' best knowledge, this is the first work to compare the bihormonal AP system and the optimal uni-hormonal AP system fairly and systematically.
The rest of this paper is organized as follows. In Section 2, the simulation models were briefly introduced. In Section 3, the IMC-PID controller was described and two switching rules were also proposed for blood glucose management. As a comparative standard, the optimal prediction/suspending therapy was presented for minimizing hypoglycemia episodes. In Section 4, the simulation results demonstrated the feasibility and safety of BG regulation using the switching control theory; statistics results can be used to further compare these therapies quantitatively; the robustness of these therapies was compared with respect to meals size variations, insulin/glucagon sensitivity variations, and measurement noises. In Section 5, to test its robustness with respect to intersubject variability, the proposed therapies were tested on ten virtual subjects. Finally, Section 6 concludes this note.

Brief Introduction on Simulation Models
The dynamic model for a virtual subject [26] consists of a glucose subsystem, an insulin subsystem, a meal subsystem, and a glucagon subsystem. The detailed descriptions for the glucose and insulin subsystems were introduced in [27,28], respectively; the meal subsystem was presented in the paper [29]. The glucagon subsystem was introduced in the paper [30].
The block diagram of the virtual subject is shown in Figure 1. 2.1. Glucose Subsystem. The glucose subsystem is described by the following four equations [27]: where is the distribution volume of the accessible compartment and 1 and 2 are the masses of glucose in the accessible and nonaccessible compartments, respectively. ( ) represents the glucose (measurable) concentration (mg/dL), 12 represents the transfer rate constant between 1 and 2 , EGP 0 is endogenous glucose production (EGP) extrapolated to the zero insulin concentration, 01 is the amount of noninsulin-dependent glucose corrected for the ambient glucose concentration, and is the rate of renal glucose clearance above the glucose threshold of 162 (mg/dL).
represents the gut absorption rate, max, is the time of maximum appearance rate of glucose in the accessible glucose compartment, and and represent the amount of carbohydrates digested and carbohydrate bioavailability, respectively.

Insulin Subsystem.
The insulin subsystem is described as follows [28]: where 1 represents mass of insulin administered as continuous infusion, 1 is the mass of insulin given as a bolus, represents the proportion of the total input flux passing through the slower, two compartment channel, 1 and 2 represent transfer rates, LD and LD are local degradation at the injection site for continuous infusion and bolus, 2 and 3 are the dose of in the nonaccessible subcutaneous compartments and the plasma compartment, respectively. ,LD represents the amount of insulin mass in which insulin degradation is equal to half of its maximal value for continuous infusion and bolus, and MAX,LD represents the saturation level describing Michaelis-Menten dynamics of insulin degradation for continuous infusion and bolus.

Meal Subsystem.
The meal subsystem is described by [29,31,32]. This model describes glucose transit through the upper small intestine and stomach. Consideṙ where sto and gut present the amounts of glucose in the stomach and intestine, respectively; sto1 and sto2 are the solid and liquid phase, respectively; ( ) presents the impulse function; presents the amount of ingested glucose; 21 and empt are the rates of grinding and gastric emptying, respectively; the abs present the rate constant of intestinal absorption.

Glucagon Subsystem.
At the heart of the model is a twocompartment representation of glucagon kinetics [30]: where ℎ is the glycogen conversion rate; max, and AG present the time of maximum appearance rate of glucose and

IMC-PID Controller.
Due to its simple structure and excellent robustness, PID controller finds widespread implementation in the industrial process control. Therefore, a great deal of effort has been directed at designing the best turning parameters for different process models. Among various PIDtuning methods, internal model control (IMC) theory has greatly succeeded because of its simplicity, good robustness, and successful practical applications [33]. For short, this kind of PID is named IMC-PID. It is well known that different patients have various dynamics, so finding suitable turning parameters is a great challenge for clinical application. One promising candidate to handle this challenge is IMC-PID, because it has only one adjustable parameter . The structure of the IMC controller is demonstrated in Figure 2, where , , and are the input, output, and load disturbances, respectively. ( ) is the controlled subject, ( ) is the subject's model, and ( ) = −1 − ( ) ( ) ( − ( ) is the minimum phase portion of ( ) and ( ) is a designed filter). In Figure 2, the part inside the dashed frame is a PID controller, , described by the following equation: From (5), one can determine the parameters of the PID controller. The order is decided by the order of the minimum phase portion ( ( )). The tradeoff between robust stability and output tracking performance can be balanced by the parameter . For smaller value of , the closedloop system will has better output tracking performance but worse robustness, vice versa. The simulation results in Section 4 demonstrate that the proposed algorithm has great robustness.

Perfect Predictive Algorithm and Insulin Pump Suspension.
According to the published studies, we have observed that the risk of "overcorrection hypoglycemia" is the greatest after a large amount of insulin. The simplest strategy to reduce hypoglycemia is to interrupt insulin delivery. It is clear that a key factor for suspending insulin infusion is having an accurate glucose predictor. Because virtual subjects were used in this study, the future glucose level can be exactly predicted using the subjects' accurate model. Hence, the optimal time for the insulin infusion suspension can be obtained such that the lowest BG value is exactly higher than 70 mg/dL.
The dosing of insulin was determined by the following equation: where 0 is the basal insulin delivery rate [35,36] and PID is designed by using IMC-PID controller. If a pending hypoglycemia event was alarmed by the perfect predictive algorithm, the insulin delivery rate will be forced to zero. Obviously, the above-mentioned algorithm can be considered the optimal uni-hormonal system.

Switched Control
System. The switched control system in this study is a dynamic system that consists of two subsystems and a switching law [37] that orchestrates switching between these two subsystems [38]. The structure of the switched controller is shown in Figure 3.
Two kinds of switching rules [39] are proposed in this study as shown in the following equations.
(1) P-type switching rule: consider where ( ) = 1 indicates insulin subsystem active mode and ( ) = 0 indicates glucagon subsystem active mode; is the blood glucose tracking error; the positive constants 1 and It shows the corresponding insulin infusion rate determined by switching controllers (P-type dual infusion therapy and PD-type dual infusion therapy, resp.). (c) It shows the corresponding glucagon infusion rate determined by switching controllers (P-type dual infusion therapy and PD-type dual infusion therapy, resp.).
(2) PD-type switching rule: condider where 1 and 2 are thresholds; > 0 is a designed parameter. Substantially, can be considered a prediction of the future tracking error , where is the prediction horizon.

In Silico Tests on Standard Virtual Subject
Based on metabolism model introduced in Section 2 and its suggested parameter settings in the literature [27][28][29][30], one virtual subject was built. For convenience, this subject was named standard subject. All simulation tests in Section 4 were finished on the standard subject. The test duration for each algorithm is 24 hours. Four kinds of closed-loop control methods were tested in each study: insulin-only therapy, prediction/suspending therapy, P-type dual infusion therapy, and PD-type dual infusion therapy. Total carbohydrate consumption consists of breakfast, lunch, and dinner (40, 60, and 85 g of carbohydrates, taking at 7:00, 12:00, and 18:00, resp.).

Nominal Case. All simulation tests consist of four cases.
There is only insulin infusion determined by = 0.42 + PID in the first case, where PID is designed by using IMC-PID controller. Figure 4(a) shows the closed-loop control performance of the insulin-only therapy, where = 0.40 was used for the tuning parameter. However, the control result is still not good enough. Hypoglycemia occurred during the closed-loop test. The second case is testing the glucose prediction and insulin suspension therapy. Through the combination of basal insulin and accurate suspending, tighter control of blood glucose can be achieved compared with the insulin-only therapy. As shown in Figure 4(a), the prediction/suspending therapy has superior capability to avoid hypoglycemia compared with the insulin-only therapy. The control algorithm can be described by  The P-type switching controller was used in the third experiment. The dosing of insulin and glucagon was determined by IMC-PID controllers, respectively. The switching rule was demonstrated by the following equations: The fourth strategy is PD-type switching therapy. The dosing of insulin and glucagon was designed by IMC-PID controllers, respectively. The switching rule was described by The excellent glycemic control results under two switching controllers (P-type and PD-type) were presented in Figure 4(a) while the designed insulin and glucagon delivery rates were shown in Figures 4(b) and 4(c), respectively. Some statistical results in Table 1 can be used to further compare the above-mentioned therapies quantitatively. The following results demonstrate that the switched therapy based on PD-type switching rule has the best glycemic control performance among these therapies.

Robustness Analysis.
The robustness of a BG management therapy is a significant factor for clinical application. Some uncertainties of external conditions may lead to severe BG fluctuation. For example, the meal size and the sensitivities of insulin, and glucagon may change dramatically.
In the sequel, the robustness analysis consists of four components: studies on meal size variations, insulin sensitivity variations, glucagon sensitivity variations, and measurement  noises. Figures 5, 6, and 7 show the BG curves under four therapies with meal size uncertainties, insulin sensitivity uncertainties, and glucagon sensitivity uncertainties, respectively, where the PD-type switching therapy has the best closed-loop control performance. The following results can help people to evaluate the robustness and safety of the proposed switching therapies. Some detailed statistical results were given in Tables 2, 3, and 4. When there are variations on meal sizes or insulin sensitivity, both insulinonly and prediction/suspending therapies are faced with increased hypoglycemia and hyperglycemia events; however, both P-type and PD-type therapies can keep all glucose concentrations within the safe range (70∼180 mg/dL). This demonstrates that the bihormonal AP system has superior robustness with respect to meal sizes and insulin sensitivity variations compared with the uni-hormonal AP system. In all six situations, PD-type therapy always has smaller BGI and SD values compared with P-type therapy. That is to say, PDtype therapy can control the glucose level more smoothly.
In clinical practice, the inaccuracy of the continuous glucose monitoring system (CGMS) enormously affects the performance of the closed-loop control system. To simulate the measurement noises, white Gaussian noises with zero mean and √ 3 standard deviation have been added on the output. The closed-loop control results with measurement noises are shown in Figure 8. One can see that all blood glucose values can be kept in the safe range even though there are measurement noises, which demonstrates that the proposed method has good robustness with respect to measurement noises. In all five situations, PD-type therapy always has smaller BGI and SD values compared with Ptype therapy. That is to say, PD-type therapy can control the glucose level more smoothly.

Method to Create Group Virtual Subjects.
In order to test whether the controller performed well for different subjects, the proposed strategy is evaluated on ten virtual subjects. Model parameters represent the high intersubject and temporal variation in insulin needs in type 1 diabetic subjects. Different combinations of free parameters in the metabolism model represent different virtual patients. These parameters are assumed to be log-normally distributed to ensure their non-negativity. Ten virtual subjects were generated using the joint distribution; that is, ten realizations of the log-transformed parameter vector were randomly extracted from the multivariate normal distribution characterizing intersubject variability. Finally, the parameters in these ten virtual subjects were obtained by using antitransformation [40]. In our opinion, these ten subjects could represent a wide range of type 1 diabetic population.  Figure 9: The blood glucose curves of ten virtual subjects under two proposed therapies: (a) P-type therapy; (b) PD-type therapy.

Simulation
Results. All simulation results were demonstrated in BG response curves and control-variability grid analysis (CVGA) plots [41]. All of these ten subjects followed a two-day scenario. The BG response curves and CVGA plot in the second day (under the proposed therapies) are given in Figures 9 and 10, respectively. From Figure 9, one can see that all subjects have no hypoglycemic event and only one tiny hyperglycemic event under two therapies. In addition, 90% of subjects are within A-zone under both P-type and PD-type therapies as shown in Figure 10. These results indicate that these proposed closed-loop control algorithms can provide satisfactory BG control, and hence it is an excellent candidate for tight blood glucose control.

Conclusions
Two novel methods based on switched control theory, Ptype and PD-type switching therapies, were proposed for  Figure 10: The control-variability grid analysis (CVGA) plot for these two proposed therapies: 90% in A-zone and 10% in B-zone for P-type therapy ( ); 90% in A-zone and 10% in B-zone for PD-type therapy (◼). the bihormonal AP system. Systematic in silico tests demonstrate that these proposed bihormonal AP systems have superior glycemic control performance than the optimal uni-hormonal AP system. Using the bihormonal AP system, all glucose concentrations can be kept within the safe range. In addition, these proposed switching therapies have excellent robustness with respect to meal size variations, insulin/glucagon sensitivity variations, measurement noises, and intersubject variability. In terms of BGI and SD, the PD-type therapy is better than the P-type one. Hence, the bihormonal AP system based on switching control theory is a promising direction for BG management.