Revitalization of cultural heritage in the digital era: A case study in Taiwan

1.   Introduction

Blood glucose needs to be maintained within a narrow range, for example, from 3.6 to 6.0 mmol/l (65 to 108 mg/dl) for the fasting people without diabetes. This narrow range differs in people and in time. A high blood glucose concentration outside the range is a symptom of diabetes, one of the prevalent diseases affecting millions of people worldwide. A peptide hormone insulin is required to interact with glucose to achieve this goal.

The interaction between blood glucose and insulin involves multiple feedback control mechanisms and is of extremely high complexity. Here we give a very basic description about it. Glucose comes from food and liver, and utilized by brain and nerve cells (insulin-independent) via the glucose transporter 3 (GLUT3) or by tissue cells such as muscle, kidney, and fat cells (insulin-dependent) via the glucose transporter 4 (GLUT4). Glucose is transported into and out of liver cells by the concentration-driven GLUT2, which is insulin-independent^[17]. In response to a low blood glucose level, the $\alpha$ cells of the pancreas produce the hormone glucagon. The glucagon initiates a series of activations of kinases, and finally leads to the activation of the glycogen phosphorylase, which catalyzes the breakdown of glycogen into glucose^[14]. In addition, the series of activations of kinases also result in the inhibition of glycogen synthase and then stop the conversion of glucose to glycogen. In response to a high blood glucose level, the $\beta$ cells of the pancreas secrete insulin^[15,35,36]. Insulin triggers a series of reactions to activate glycogen synthase, which catalyzes the conversion of glucose into glycogen^[16,21]. Insulin also initiates a series of activations of kinases in tissue cells to lead to the redistribution of the glucose transporter 4 (GLUT4) from intracellular storage sites to the plasma membrane^[13]. Once at the cell surface^[8,32], GLUT4 transports glucose into the muscle or fat cells.

Mathematical models for the blood glucose regulation system are required for integrating a glucose monitoring system into insulin pump technology to form a closed-loop insulin delivery system, the so-called "artificial pancreas" (see Hovorka^[18], Panteleon et al.^[33], Steil et al.^[38]). To make this artificial pancreas as close as possible to the natural pancreas, many mathematical models about the interaction between blood glucose and insulin have been proposed. In many of these models, the compartmental modeling approach has been used. In the 1960's, Ackerman et al.^[1,2] proposed a simplified compartmental model of two linear ordinary differential equations for the glucose-tolerance test. This model lumped the large number of kinetic parameters into a much smaller number which can, at least in part, characterize the human glucose regulatory system. In the 1970's, Bergman et al.^[4,5,6,7] proposed a minimal compartmental model of nonlinear ordinary differential equations. It was assumed that the blood glucose and hormone insulin are contained in two different compartments and interact with each other. The plasma insulin is secreted from the endocrine system under the stimulation of the blood glucose and enters a "remote compartment", where it is active in accelerating glucose utilization. Then the dynamics of the blood glucose, plasma insulin, and the insulin in the remote compartment was modeled by a system of three nonlinear ordinary differential equations. In the 1990's, to determine whether the oscillations of insulin and blood glucose could result from the feedback loops between insulin and glucose, Sturis et al.^[39,40] developed a parsimonious model including the major mechanisms involved in glucose regulation. This model comprises two major negative feedback loops describing the effects of insulin on glucose utilization and glucose production, respectively, and both loops include the stimulatory effect of glucose on insulin secretion. Then this model consists of six nonlinear ordinary differential equations, three of which are used to model a time delay for the effect of insulin on glucose production and a sluggish effect of insulin on glucose utilization. On the basis of Sturis et al.'s model, Li et al.^[23] introduced two explicit time delays and proposed a more robust alternative model for better understanding the glucose-insulin endocrine metabolic regulatory system and the ultradian insulin secretory oscillations for the cases of continuous enteral nutrition and constant glucose infusion. More complicated compartmental models proposed by Hovorka et al.^[18], Liu et al.^[24,25,27], Man et al.^[29,30,31], and Sorensen^[37] have considered labeling more compartments for better understanding the behavior of different parts of the body. To assess the liver ability to extract insulin, Toffolo et al.^[41,42], proposed a minimal model of insulin delivery and kinetics by combining the classical C-peptide minimal model. In the 2010's, Vahidi et al.^[45,46] developed a model for type Ⅱ diabetes, based on a detailed compartmental physiological model proposed by Sorensen^[37] for a healthy human body. There have been a huge set of references on the modeling research and it is hopeless to mention all of them.

In all these models, no parameter uncertainties were introduced. The parameter uncertainties here means that a parameter may vary in a neighborhood of its nominal value. This can happen in practice. For instance, an estimation error may be introduced when a model parameter is estimated. In addition, parameter values may vary slightly from person to person. Therefore, it is of practical interest to introduce parameter uncertainties into a model. As far as the robust tracking is concerned, the robustness against the parameter uncertainties means that a feedback and feedforward controller should be designed in a way such that it can tolerate the parameter variation in a small neighborhood of their nominal values. That is, even though the parameters have a small change, the controller can still achieve its tracking goal. Such a robust problem may be solved by only a robust control design method and other methods such as the sensitivity analysis, bootstrapping to quantify parameter ranges, and stochastic modeling might not be relevant.

As a first attempt in this direction, we introduce parameter uncertainties into the model proposed by Bergman et al.^[4,5]. Although this model is simple and consists of only three differential equations, one equation for the glucose dynamics, one equation for the effect of the remote insulin on glucose, and the other equation for the hormone insulin dynamics, it can simulate well the major biological dynamics about the interaction between glucose and insulin, as evidenced in the numerous literature (see. e.g., ^[4,5,6,7,11]). It has been widely used and further updated later (see, e.g., ^{[10,18,20,23,24,39,40,41]}).

Furthermore, we also include a time-dependent exogenous glucose input $J$ from food through the intestine in this model and assume that the desired glucose reference $r$ to be tracked is time-dependent, as described in Figure 1. Then we design a robust feedback (using an actual glucose concentration $g$) and feedforward (using the exogenous glucose input $J$ and the desired glucose reference $r$) controller $u(g, r, J)$ such that the glucose concentration asymptotically tracks the time-varying glucose reference for any small uncertainties. Numerical simulations conducted with published experimental blood glucose data show that our controller is effective in achieving the blood glucose tracking and agree with our theoretical results.

2.   Formulation of a model with parameter uncertainties

In the 1970's, Bergman et al.^[4,5] proposed a mathematical model for the blood glucose regulation system. To establish such a model, they made a simplified assumption. The blood glucose is produced from food through the intestine. The blood glucose and hormone insulin are contained in two different compartments and interact with each other. The blood glucose can be converted into the liver glycogen and used by the tissue metabolism with or without the help of insulin. The plasma insulin is secreted from the endocrine system under the stimulation of the blood glucose and enters a "remote compartment", where it is active in accelerating glucose utilization. It is used up by the insulin metabolism. Under this assumption, the dynamics of the blood glucose and insulin was modeled by the following system of nonlinear differential equations ^[5]:

In the above equations, $g$ and $h$ denote concentrations of blood glucose and plasma insulin, respectively, $h_a$ is the effect of the remote insulin on glucose, $J$ is a rate of the exogenous glucose input from the intestine, $u$ is a rate of insulin secreted from the endocrine system or infused externally, and $m_1, m_2, m_3, m_4$ are positive rate constants. The biological meaning of these rate constants are explained in the Table 1. Because the glucose needs to be monitored and regulated, we introduce an output equation:

To maintain the blood glucose at a normal level of about 6 mmol/l, a basal exogenous glucose input $G_b$ is required^[24,47]. Thus the blood glucose input $J$ could be split into two parts:

where $G_d$ can be treated as a glucose disturbance.

In what follows, $(\bar{g}, \bar{h}_a, \bar{h}, \bar{u})$ denotes the equilibrium of the system (2.1)-(2.3) without the glucose disturbance. That is, they are solutions of the following steady state system:

In practice, an estimation error may be introduced when a model parameter is estimated. In addition, the values of the parameters $m_1, m_2, m_3, m_4$ may vary slightly from person to person. Therefore, it is reasonable to introduce parameter uncertainties into the system (2.1)-(2.3):

In the above equations, $m_1, m_2, m_3$ and $m_4$ represent the nominal part of the parameters and $w_1, w_2, w_3, w_4$ represent the uncertain part. Since the equilibrium $(\bar{g}, \bar{h}_a, \bar{h}, \bar{u})$ would not depend on the uncertainties, the terms $w_1 \bar{g}, w_2 \bar{h}_a, w_3 \bar{h}$ and $w_4\bar{h}$ are included such that this is guaranteed.

Because the blood glucose needs to be maintained within a narrow range from 3.6 to 6.0 mmol/l (65 to 108 mg/dl), we need to develop an insulin infusion rate $u$ in response to a glucose surge to achieve this goal. This problem can be mathematically formulated as follows. Given a glucose reference $r(t)$ (a constant 6 mmol/l or fluctuating around 6 mmol/l), we want to design a controller $u = u(g, J, r)$ such that for all sufficiently small uncertainties $w_1, w_2, w_3, w_4$

3.   A static feedback and feedforward controller

To solve this control problem, we need to express $G_d$ and $r$ by using solutions of linear differential systems. The glucose produced from food via the intestine is not exactly known and might be random. This glucose disturbance could be of periodic nature and then can be approximated by a Fourier polynomial

The number $1$ and functions cosine's and sine's can be generated by the linear differential system

where

Then $G_d$ can be expressed as

where $\mathbf{C}_d = [a_0, a_1, b_1, a_2, b_2, \cdots, a_n, b_n]$ is a constant vector.

Since the glucose reference $r$ is of periodic nature, too, it can be generated in the same way as follows:

where $\bar{g}$ is the glucose equilibrium, $\mathbf{C}_r$ is a constant vector and $\mathbf{v}_r$ is governed by

with $\mathbf{A}_{r}$ being similar to $\mathbf{A}_{d}$:

Let

Then glucose reference and disturbance can be lumped together as follows:

We introduce the tracking error:

Combining the exosystem (3.2) with the model (2.9)-(2.11), we obtain a composite system

When the disturbance $G_d = 0$ and $r = \bar{g}$, we can readily show that $g, h_a, h$ and $u$ converge to their constant equilibrium, respectively. When $G_d$ and $r$ are time-dependent, we can expect that $g, h_a, h$ and $u$ would converge to their time-varying "equilibrium" that would depend on $\mathbf{v}$. In fact, such time-varying "equilibrium" called a stable center manifold is guaranteed by the center manifold theory (see, e.g., ^[9]). Thus we introduce the following $\mathbf{v}$-dependent variable transform:

to translate the time-varying "equilibrium" to the origin. Substituting the above transformation into the equations (3.3)-(3.7), we obtain

where $\mathbf{w} = [w_1, w_2, w_3, w_4]$ and

In order to eliminate the terms that do not contain either $\hat{g}$ or $\hat{h}$, we set

Notice that the equations (3.17)-(3.19) are the result of the center manifold theory (see^[9]). In this way, the tracking problem (3.3)-(3.7) is converted to the following stabilization problem:

The regulator equations (3.17)-(3.20) can be solved explicitly as follows:

It is easy to see that the equilibrium $(0, 0, 0)$ of the system (3.21)-(3.24) at $\mathbf{v} = \mathbf{0}$ can be stabilized exponentially and locally by a static output feedback

Therefore, we have constructed a static output feedback

to regulate the blood glucose to its desired reference $\bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}$. Incorporating this feedback controller into the system (3.3)-(3.7), we obtain a closed-loop system:

We summarize this result in the following theorem.

Theorem 3.1. Suppose that $k_g$ is a constant to exponentially stabilize the equilibrium $(0, 0, 0)$ of the system (3.21)-(3.24). Then the closed-loop nonlinear system (3.32)-(3.36) asymptotically tracks the blood glucose reference $\bar{g} +[\mathbf{C}_r, \; \mathbf{0}] \mathbf{v}(t)$, that is,

In numerical simulations below, we use the software Matlab and the experimental data from Sturis et al.'s ^[39]. We use the software Engauge Digitizer 4.1 to read the data.

To conduct a numerical simulation, we use the experimental data from Sturis et al.'s ^[39] to determine $\mathbf{C}_d$, $\mathbf{C}_r$, and $\mathbf{A}_e$. These data were obtained from normal men without a history of diabetes and the detailed experimental methods and procedures are referred to^[39]. Since the data contain 24 hours profiles of glucose and plasma insulin, they are appropriate for testing whether the closed-loop nonlinear system (3.32)-(3.36) asymptotically tracks the blood glucose reference. The data were read from ^[39] by using the software Engauge Digitizer 4.1. According to the second panel of Figure 2 of Sturis et al.'s ^[39], the exogenous glucose infusion shown in Figure 2 is given by

Because the glucose concentration unit is mM, we need to covert the unit of cc per hour of the infused solution to the unit of mM per minute of the infusion rate inside the body. According to Sturis et al.'s paper^[39], the infused glucose solution is 20% solution. This could mean that one cc solution contains 0.2 grams of glucose. The blood volume of a person with the weight of 70 kg is about 6 liter. Using the conversion 1 mM = 18 mg/dl, we obtain

In this simulation, we find that $G_b = 0.3218$ mM/min. Thus we obtain

and

Fitting the experimental glucose data shown in the second panel of Figure 2 of Sturis et al.'s ^[39] into the Fourier polynomial:

we obtain

The unit of the coefficients $a_0$, $a_i$ and $b_i$ is mM, where $i = 1, 2, \cdots, 17$. The glucose equilibrium is set to $\bar{g} = 6$ mM. Thus we obtain

and

Figure 3 shows that the data can be fitted into this Fourier polynomial very well.

The estimation of the parameters $m_1, m_2, m_3$, and $m_4$ is challenging. Because the expression of the controller $u$ is too big to be expressed in Matlab, we failed to use the Matlab function "sbioparamestim" to estimate them. Thus, starting with the values from Bergman et al.'s paper^[5]: $m_1 = 0.0296$ (/min), $m_2 = 0.0186$ (/min), $m_3 = 0.0000065$ (/($\mu$U/ml)/min$^2$), $m_4 = 0.238$ (/min), we manually adjust them to achieve as a best result as we can. These estimated values are presented in Table 1. Compared with Bergman et al.'s original values, only $m_1$ has a large change. The uncertainties $w_1, w_2, w_3,$ and $w_4$ are taken to be 10% of their nominal values. The proportional constant $k_g$ is taken to be $0.01$ (/min).

The model is solved numerically with Matlab. We first run the model 1500 minutes to achieve a stable state and then set that time as the initial time. Figure 4 shows that the blood glucose $g$ is asymptotically tracking the desired time-varying reference determined by the experimental data from Sturis et al.'s ^[39]. Although the simulated insulin profile does not exactly match with the experimental data, the simulated and experimental pulses occur simultaneously. The dynamics of the controller $u$ and insulin effect $h_a$ follows the dynamics of the glucose infusion.

Because the static output feedback (3.31) contains the uncertainties $\mathbf{w}$, it is not robust. If the uncertainties $\mathbf{w}$ are set to 0 in the static output feedback controller (3.31), Figure 5 shows that the glucose $g$ fails to track its reference. In fact, it has been shown that there are no static state feedback controls that can solve the robust tracking problem (see^[19]).

4.   A robust dynamical feedback and feedforward controller

We now design a robust dynamical feedback and feedforward controller via linearization. The linearization of the system (3.3)-(3.7) at its equilibrium is given by

where

To construct a robust control law, we introduce the following dynamic compensator:

where the matrix pair $(\mathbf{G}_1, \mathbf{G}_2)$ is determined below. Since the state variables $\tilde{h}_a$ and $\tilde{h}$ are not available for feedback, we also need to use the Luenburger observer:

where the vector $\mathbf{L}$ is to be designed and

We then introduce the following dynamic output feedback and feedforward controller:

where the vectors $\mathbf{K}_o$ and $\mathbf{K}_z$ are to be designed.

A minimal polynomial of a matrix $\mathbf{A}$ is a polynomial $P(\lambda)$ with the least degree such that $P(\mathbf{A}) = \mathbf{0}$. If $\omega_r = \omega_d$, then the minimal polynomial of $\mathbf{A}_{e}$ is given by

where $l = \max(m, n)$. If $\omega_r\neq\omega_d$, then the minimal polynomial of $\mathbf{A}_{e}$ is given by

where $l = m+n$. Then $\mathbf{G}_1$ is defined to be a companion matrix of the minimal polynomial and $\mathbf{G}_2$ is set to be a matrix such that the pair $(\mathbf{G}_1, \mathbf{G}_2)$ is controllable as follows:

In control theory, such a matrix pair $(\mathbf{G}_1, \mathbf{G}_2)$ is called to incorporate a 1-copy internal model of the matrix $\mathbf{A}_{e}$. In the following proof, we can see the reason why the pair $(\mathbf{G}_1, \mathbf{G}_2)$ is set up this way.

To convert this tracking problem to a stabilization problem, we introduce the coordinate transformation

Substituting this transformation (4.12) and the feedback (4.10) into the composite system (4.1)-(4.5) and the equations (4.7) and (4.8), we obtain

In order to eliminate $\mathbf{v}$ from the above system, we set

and then

If the regulator equations (4.13)-(4.18) have a solution, then the tracking problem is converted to the problem of stabilization of the system (4.19)-(4.24).

We first show that the equilibrium $(0, 0, 0, \mathbf{0}, \mathbf{0})$ of the system (4.19)-(4.24) can be exponentially stabilized at $\mathbf{w} = \mathbf{0}$ and then robustly stabilized for sufficiently small $\mathbf{w}$. For this, we set $\hat{\mathbf{x}} = [\hat{g}, \hat{h}_a, \hat{h}]^T$ and rewrite the system (4.19)-(4.24) at $\mathbf{w} = \mathbf{0}$ as follows:

Setting $\hat{\mathbf{q}} = \hat{\mathbf{x}}-\hat{\mathbf{o}}$, we derive from the above system that

and then

It is easy to see that $(\mathbf{A}, \mathbf{C})$ is observable and the pair

is stabilizable (see, e.g., ^[19]). So there exist constant vectors $\mathbf{K}_o, \mathbf{K}_z, \mathbf{L}$ such that all eigenvalues of the coefficient matrix of the system

have negative real parts. That is, the equilibrium $(0, 0, 0, \mathbf{0}, \mathbf{0})$ of the system (4.19)-(4.24) can be exponentially stabilized at $\mathbf{w} = \mathbf{0}$.

So it remains to show that the regulator equations (4.13)-(4.18) have a solution. To this end, we set $\mathbf{X} = [\mathbf{G}, \mathbf{H}_a, \mathbf{H}]^T$ and rewrite the equations (4.13)-(4.17) as one matrix equation:

where

Since $\mathbf{Q}$ and $\mathbf{\Omega}$ are equivalent, all eigenvalues of $\mathbf{Q}$ also have negative real parts. Because all eigenvalues of $\mathbf{A}_e$ have zero real parts, the above matrix equation has a unique solution (see^[19]). To show the solution also satisfies the equation (4.18), we let $\mathbf{Z} = [\mathbf{Z}_{1}, \mathbf{Z}_{2}, \cdots, ,\mathbf{Z}_{2l+1}]^T$. It then follows from the equation (4.16) that

Here we have used the property that $p(\lambda)$ is the minimal polynomial of $\mathbf{A}_e$, that is, $p(\mathbf{A}_e) = \mathbf{0}$. The above proof of the equation (4.18) gives the reason why the pair $(\mathbf{G}_1, \mathbf{G}_2)$ is constructed as in (4.11).

Using the translation (4.6), we obtain a dynamic output feedback and feedforward controller for the original nonlinear system (3.3)-(3.7):

For the convenience of reference, we collect all relevant equations together to form a closed-loop system:

Theorem 4.1. Suppose that $\mathbf{K}_o, \mathbf{K}_z, \mathbf{L}$ are constant vectors such that the matrix

is exponentially stable. Then the blood glucose of the closed-loop nonlinear system (4.43)-(4.49) asymptotically tracks its reference $\bar{g} + [\mathbf{C}_r, \; \mathbf{0}] \mathbf{v}(t)$ for any sufficiently small uncertainty $\mathbf{w}$, that is,

We use the same data as used in the above section to conduct a numerical simulation for this robust approximate controller (4.40). In this simulation, we first used an internal model 1-copy $\mathbf{G}_1$ of the matrix

Unfortunately, Matlab failed to solve the closed-loop control system (4.43)-(4.49) with this high dimensional $\mathbf{G}_1$. Thus we have to reduce the dimension to obtain an internal model 1-copy $\mathbf{G}_1$ of the matrix

In this case, the matrices $\mathbf{G}_1$ and $\mathbf{G}_2$ are given by

The vectors $\mathbf{K}_o$ and $\mathbf{K}_z$ are determined in the way such that the eigenvalues of the matrix $\mathbf{\Phi}- \mathbf{\Psi}\mathbf{K}$ are: $-0.0500, -0.0909, -0.1318, -0.1727, -0.2136, -0.2545, -0.2955, -0.3364,$ $-0.3773, -0.4182, -0.4591, -0.5000,$ where $\mathbf{K} = [\mathbf{K}_o, \mathbf{K}_z]$. Using the algorithm for computing the feedback gain matrix from Section 3.4 of the book^[26], we obtain

The vector $\mathbf{L}$ is determined in the way such that the eigenvalues of the matrix $\mathbf{A}- \mathbf{L}\mathbf{C}$ are: $-0.01, -0.02, -0.03$. Using the Matlab function place, we obtain

According to our simulation studies, the eigenvalues should be between $-0.5$ and $-0.01$ to achieve a good performance of the controller. All other parameter values are the same as in the previous simulations.

Comparing Figure 6 with Figure 4, we found it interesting that the dynamics of all state variables under the robust controller (4.40) is identical to the one under the static controller (3.31).

The static controller (3.31) and the robust dynamic controller (4.40) have their own advantages and disadvantages. The static controller (3.31) is easy to implement, but is not flexible, no parameters being used to adjust to improve its performance. Since it is not robust, a small change in the model parameters results in a deviation of blood glucose from its target reference. On the other hand, the robust dynamical controller (4.40) can tolerate certain small changes in the model parameters and its performance can be improved by adjusting certain controller parameters such as $\mathbf{K}_o, \mathbf{K}_z$ and $\mathbf{L}$. However, it is more difficult to implement because the dimension of the closed-loop control system is greatly increased through the dynamic compensator (4.7). As discussed above, Matlab failed to solve the closed-loop control system numerically if the original $\mathbf{G}_1$ for the matrix $\mathbf{A}_e$ is used. Furthermore, its performance critically depends on the gain vectors $\mathbf{K}_o, \mathbf{K}_z$ and $\mathbf{L}$. It is a great challenging problem to determine them to achieve the best performance of the controller in a way such that the simulated both glucose and insulin dynamics agrees with the real experimental data. Since the controller performance issue itself is an independent research area^[12,22,28], this problem is beyond the scope of this work.

5.   Discussion and conclusion

In this work, we have solved a robust glucose tracking problem. Parameter uncertainties were introduced into the mathematical model proposed by Bergman et al.^[5]. Using an actual blood glucose level as feedback and an exogenous glucose input and a desired glucose reference as feedforward, we have designed a robust feedback and feedforward controller such that the blood glucose tracks the desired time-varying glucose reference for any small uncertainties. Numerical simulations conducted with published experimental blood glucose data show that our controller is effective in achieving the blood glucose tracking and agree with our theoretical results.

Limitations exist in this work. Although our robust feedback and feedforward controllers might have a potential application in the insulin pump technology, more comprehensive tests and improvements with experimental and clinic data are needed to be done before they might be implemented in an insulin pump. The interaction between glucose and insulin is complex. Our current simple model missed considering many other important control mechanisms such as the glucose regulation by the hormone glucagon and the organ liver. Therefore we may need to consider more complex mathematical models, such as the model proposed by Sturis et al.^[40] and the model proposed by Liu et al.^[24], and introduce parameter uncertainties into them. In all these cases, the regulator equations are more complex and it might be impossible to solve them analytically. Then a different approach such as used in^[3] is needed to design a controller.

Acknowledgments

The author heartily appreciates the reviewers' constructive comments, which greatly enhance the quality of the paper to meet the journal standard.

Conflict of interest

The author declares no conflicts of interest in this paper.