January 3, 2022

By Hannah Le & Lucy Ma under the guidance of Dr. Paul Yoo

The typical routine for diabetes management includes injecting insulin to a general area of the body, anywhere from one to four times per day, which can be often painful and uncomfortable for patients. In fact, a study in 2016 reveals that barriers to insulin therapy adherence include injections being time consuming and lacking sufficient injection instructions. We would like to propose a minimally invasive close-loop control method for continuously delivering insulin using a PID controller, modeled by Simulink. This paper proposes a Proportional Integral Derivative (PID) model that controls and maintains the blood glucose level through controlling insulin delivery. This paper used the Bergman Minimal model to predict how blood glucose concentration fluctuates in response to insulin injection in Type I diabetes patients. The initial tuning parameters were found using the First-Order Linear Dead Time Model (FODT), and the IMC correlations were used to manually tune and optimize the PID controller and evaluate it against the design criteria based on a physiological response to insulin injections in Type 1 diabetes patients. Using the transfer function derived from system identification and the PID transfer function, we were able to calculate our open loop gain and determine system stability through a critical analysis of the Root Locus Plot and the Nyquist Plot. The presented model can be further improved by decreasing the settling time and reducing percent overshoot to reduce the risk of hypoglycemia due to a high insulin input. Finally, with insulin injection automated based on blood glucose level input, the paper proposes further development to give the patient more control over the device by having a front-end user interface that allows the patients to approve the insulin delivery. Physicians can also take advantage of this model to receive real-time upon PID automatic detection.

In the human body, the pancreas is responsible for sensing blood glucose concentration, calculating the amount of insulin needed, and then delivering the necessary amount to restore normal blood glucose levels. In this study, we propose to model a device that acts as an “artificial pancreas” consisting of a glucose sensor, insulin pump, and a control algorithm used to regulate the pump to deliver insulin in order to maintain a normal blood glucose level. This device is preferred because it avoids the need for finger pricks and blood draw, numerous insulin injections, and reduces the risk of hypoglycemia/hyperglycemia. We model the closed-loop control process that takes place in blood glucose regulation by insulin. This project aims to demonstrate the superiority of automatic control compared to manual control, generates dynamic data for modeling with first principles, and demonstrates the tunability of the proportional integral derivative (PID) controller. We hypothesize that the PID controller will be able to minimize percentage overshoot, reduce settling time and peak time, and exhibit stability through the critical analysis of the Root Locus, Bode, and Nyquist plots.

**A1** shows the diagram of the closed-loop PID control system for insulin control to regulate blood glucose levels. The insulin input to the patient is dependent on the difference between the value from the feedback loop and the desired blood glucose concentration from a reference input. This generates an error signal that is fed to the PID to determine the amount of insulin needed to be delivered to the patient to regulate the blood glucose concentration. The meal disturbances will also be an input to the model, to show the relationship between insulin and glucose concentration given the disturbances in a patient. This dynamic model will be represented using nonlinear ordinary differential equations as shown in **A2.** Using system identification, the transfer function can be derived.

The PID controller represents an insulin pump that continuously delivers insulin into a patient’s body. The controller calculates the insulin delivery based on:

- Proportional term (P), which adjusts the insulin delivery depending on the current blood glucose concentration measured
- Integral term (I), which adjusts insulin delivery depending on the area under the curve between the set point and the measured blood glucose concentration
- Derivative term (D), which delivers insulin depending on the rate of change of blood glucose over time

The set point for blood glucose concentration is a range based on a healthy individual, between 65 and 105 mg/dL with a target of 80 mg/dL, as in A3 shown by the dotted black lines of the second graph. Values outside of the safe range are not desirable as it can lead to hyperglycemia and hypoglycemia. The metabolic reactions in the body to insulin injections is shown in A4. This shows the insulin pharmacodynamics of the various commercial insulin options, and the design criteria will be derived based on the onset, peak, and duration of effect. In addition, the PID model should avoid reaching hyperglycemia levels (>105mg/dl) or hypoglycemia (< 65mg/dl). This gives the model a “safe range” to tune the parameters to. The rise time should be short, the steady-error should be eliminated, the overshoot minimized, the settling time < 5 hours, and peak time < 2 hours. The quantitative values are based on values found in A4.

**Linearizing the first-principle, non-linear diabetes dynamic system with a transfer function**

To predict how blood glucose concentration fluctuates in response to insulin injection in Type I diabetes patients, we used the Bergman Minimal Model, which takes into two inputs: insulin infusion rate (mU/min) and carbohydrate intake (mmol/L-min). Such a model is appropriate for our system because we want to investigate how plasma glucose concentration changes over 24 hours with meal disturbances to various insulin inputs. One of the main challenges of using the Bergman model is it is non-linear, as seen in Equation X, which makes it difficult for us to calculate the open loop gain which is defined as the multiplication of the system gain and the control gain in the open loop. Hence, we utilized the MATLAB System Identification Toolbox to linearize the Bergman model.

System Identification Toolbox is a MATLAB app that can construct a transfer function from measured input-output data [**A5**]. At a high level, MATLAB aims to minimize the mean squared error between the data generated by MATLAB’s transfer function and Bergman’s output data using the Instrument Variable Approach (IV). It also allows one to customize the number of zeroes and poles to evaluate which transfer function performs best. Leveraging these functionalities, we followed the subsequent steps to find the transfer function of best fit:

- As glucose response diabetes has 6 states, we have 6 degrees of freedom. Hence, we can start with 1 zero and 6 poles.
- Based on the transfer function’s performance with the test dataset, we modified the number of zeroes and poles to generate a total of 35 transfer functions.
- To minimize overfitting, we apply our 35 transfer functions to 5 different validation datasets and were able to arrive at a transfer function that performs consistently well (>43% accuracy) for all data scenarios
**[A6]**. Unfortunately, we were not able to obtain any transfer function with a higher accuracy level, which was expected because the process of manual tuning is time-consuming and the Bergman nonlinear system was a sixth order system. However, we observed that the transfer function often deviated from the Bergman’s model by only 5 mg/dL of plasma glucose concentration, which was shown to be non-fatal to the human body.

**Finding initial tuning parameters using the First-Order Linear Dead Time Model (FODT)**

A first-order with dead time is a simple approximation of the dynamic response of a variable (in our case blood glucose level) to an external influence. Mathematically, the model can be stated as a linear first order equation pdy(t)dt=-y(t) + Kpu(t-p) where Kp, p, and p are process gain, process time constant, and process dead time respectively. In order to derive these three parameters, we first obtained a step glucose response to step insulin input as seen in A7 using the Bergman model. Process gain was then calculated by dividing y in step response by u in step input, giving us -36.2. We were also able to find p, the apparent dead time which was 0.1 hours by reading the step response graph. Finally, following step 4, 5, and 6 from A8 we found the process time constant p to be 1.8 hours.

**Using IMC Correlations and Manual Tuning to Optimize PID**

Using the process gain, dead time and time constant, we are able to achieve the initial tuning parameters using the IMC tuning guidelines, as seen in **A9.** Hence, our proportional, integral, and derivative gains were -0.0276, 1.8, and 0, respectively. When we tested this transfer function for insulin infusion rates input with meal disturbances, we observed that there was a steady state error of around 3 mg/dL from our target 80 mg/dL, and an overshoot of 10 mg/dL above 105 mg/dL, which was the upper limit of our safe range for blood glucose level, potentially leading to hyperglycemia. Since the derivative controller measures the rate of change of error overtime and informs us whether the system is rapidly increasing towards a set point, we expected that increasing the derivative gain beyond 0 would help. Further, since having large integral terms might lead to higher percent overshoot, we might want to slightly decrease this value while still keeping it there to decrease steady state error. Using the same proportional gain, a derivative gain of 1, and an integral gain of 0.5, we were able to achieve a response to insulin infusion rate as seen in **A10**. Here, the fluctuations of blood glucose concentration is within the safe range of 65 to 105 mg/dL. Similar to how the body responds to needle insulin injections, we achieved a rise time of 30 minutes, a peak time of 1 hour, a settling time of 4 hours and a percentage overshoot of 16.6%.

**Critical Evaluation of System Stability Using Root Locus, Bode, and Nyquist Plots**

Using the transfer function derived from system identification and the PID transfer function, we were able to calculate our open loop gain, which was seen in **A11**. From the pole-zero map in **A12** and root locus plot in **A13**, we observed that all 6 poles were on the left hand side of the imaginary axis, meaning our system is stable. Since we have three zeros on the right hand side of the imaginary axis and the Nyquist plot in **A14** to enclose -1 in a circle clockwise direction three times, we ended up with P=Z-N=0 on the RHS, which once again verifies the stability of our system. Using the Bode plot in **A15**, we can determine the frequency response of the transfer functions, and quickly communicate how the magnitude and phase of the transfer function fluctuates at different frequencies. This plot can be used to compare a PI and a PID controller, where the derivative term can be expressed as the transfer function of a phase-lead filter, which can be designed and transformed into a derivative plus filter combination. To see the effect of the phase-lead filter, the open-loop frequency response can be analyzed to see whether the filter adds phase to the loop, and increases the overall phase margin of the closed-loop [**A16**]. Another application of the Bode plot is to find the Bode’s integrals for slope adjustment in PID controller design [**A17**], which shows the relationship between the phase and the amplitude of minimum phase stable systems. This is useful because if we wanted to improve the controller performance by adjusting the phase margin and the slope of the Nyquist curve at the crossover frequency, the derivatives can be approximated by the Bode’s integrals. To obtain a desired phase margin **Φd** at the crossover frequency wc, we solve the equations shown in **A17**.

This model demonstrated the usage of a PID controller to control the delivery of insulin based on blood glucose concentration levels. Using the IMC tuning guidelines, the initial tuning parameters were achieved. Our proportional, integral, and derivative gains were -0.0276, 1.8, and 0, respectively. There was a steady state error of around 3 mg/dL from our target 80 mg/dL, and an overshoot of 10 mg/dL above 105 mg/dL. Since this exceeded the limit for being at risk for hyperglycemia, the parameters need to be adjusted to minimize overshoot. Using the same proportional gain, a derivative gain of 1, and an integral gain of 0.5, where the blood glucose concentration is within the safe range of 65 to 105 mg/dL as specified by our design criteria. To demonstrate that the model represents the physiological response well, we found that the model achieved a rise time of 30 minutes, a peak time of 1 hour, a settling time of 4 hours and a percentage overshoot of 16.6%, all of which meet the quantitative criteria. Through a critical analysis of the root locus and nyquist plots, we found that the system was stable. Other linearization techniques can be explored and find a better transfer function using data from first principles dynamic system for one that fits best with system identification to reduce overfitting and improve generalization to a variety of patient glucose input data.

From analyzing the output models of how the blood glucose concentration changes with different input models and parameters, it is important to tune the parameters further to minimize percent overshoot to avoid a state of hyperglycemia in patients. A large settling time and overshoot can cause the glucose concentration to overdose and lead to hypoglycemia because of high insulin input. We should also explore other linearization techniques and find a better transfer function that generalizes well to a variety of input data, using data from first principles dynamic systems for one that fits best with system identification. This avoids the problem of using a transfer function that fits well with one scenario and does not fit for others.

We can further develop this project to show an end-to-end demonstration of a remote PID controlled insulin delivery device. The device would have two functional units: a glucose sensing unit and a therapeutic unit. Lee *et. al* suggested the device uses a layer of fluoropolymer Nafion, a material effective for absorbing sweat and carrying it towards the device’s glucose sensor, which is built on modified graphene. The glucose sensor includes enzyme glucose oxidase, which can react with sugar to produce hydro peroxide. Through a series of electrochemical events, hydroperoxide extracts current from the doped graphene, creating an electrical signal proportional to the amount of glucose present. The patch can then send the signals to a wireless device that monitors them to analyze the glucose signal to determine the appropriate next steps. **[A18]**

If hyperglycemia indeed takes place, the external device sends signals wirelessly back to the therapeutic unit. This unit contains hydrogel microneedles that can go through the epidermal layer of the skin towards the bloodstream. Microneedles are known to be flexible and harmless, enabling them to sit on the skit for an extended period of time. The needles can contain diabetes drugs, such as metformin, and a dissolvable polymer such as polyvinyl pyrrolidone. Sitting on top of the microneedles is a gold and graphene mesh. They serve as a heater to melt the microneedle coating. Once an hyperglycemia event occurs, the mesh heats up tridecanoic acid, causing it to melt. Polyvinyl pyrrolidone is then dissolved and the drug payload is released. **[A18]**

The blood glucose output of the model is connected real time to a patient’s app on the phone, which gives the patient more control over the device by having a front-end user interface that allows the patients to approve the insulin delivery, upon PID automatic detection. Next, insulin delivery and blood glucose levels detected are documented and sent to physicians for patient record via real time signals processing. Given the challenges posed by COVID-19 for physicians to check on patients, this medical device is an innovative solution.

**A1**

**A2**

**A3. Blood Glucose Level Response With PID Controller and Meal Disturbances To Insulin Infusion Over 24 hours**

**A4**

Donner T, Sarkar S. Insulin – Pharmacology, Therapeutic Regimens, and Principles of Intensive Insulin Therapy. [Updated 2019 Feb 23]. In: Feingold KR, Anawalt B, Boyce A, et al., editors. Endotext [Internet]. South Dartmouth (MA): MDText.com, Inc.; 2000-. Available from: https://www.ncbi.nlm.nih.gov/books/NBK278938/

**A5**

System Identification Toolbox. (2020). Retrieved December 10, 2020, from https://www.mathworks.com/products/sysid.html

**A6**

**Simulation of TF34’s Output In Comparison with Non-linear Diabetes Model’s Output**

**Five Different Insulin Inputs to Generate the TF and Diabetes Model’s Output**

**A7**

**Step Blood Glucose Level Response To Step Insulin Input**

**A8**

Kokate, R. D., Waghmare, L. M., & Deshmukh, S. D. (2010, October). Review of tuning methods of DMC and performance evaluation with PID algorithms on a FOPDT model. In 2010 International Conference on Advances in Recent Technologies in Communication and Computing (pp. 71-75). IEEE.

**A9**

Sánchez, H. S., Padula, F., Visioli, A., & Vilanova, R. (2017). Tuning rules for robust FOPID controllers based on multi-objective optimization with FOPDT models. ISA transactions, 66, 344-361.

**A10**

**Blood Glucose Level Response With PID Controller and Meal Disturbances To Insulin Infusion Over 24 hours**

**A11**

**A12**

**A13**

**A14**

**A15**

**A16**

Kikuuwe, R., Kanaoka, K., Kumon, T., & Yamamoto, M. (2015). Phase-lead stabilization of force-projecting master-slave systems with a new sliding mode filter. IEEE Transactions on Control Systems Technology, 23(6), 2182-2194.

**A17**

Karimi, A., Garcia, D., & Longchamp, R. (2002, May). PID controller design using Bode's integrals. In Proceedings of the 2002 American Control Conference (IEEE Cat. No. CH37301) (Vol. 6, pp. 5007-5012). IEEE.

**A18**

Lee, H., Choi, T., Lee, Y. et al. A graphene-based electrochemical device with thermoresponsive microneedles for diabetes monitoring and therapy. Nature Nanotech 11, 566–572 (2016). https://doi.org/10.1038/nnano.2016.38