Mathematical modelling of cardiac pulse wave reflections due to arterial irregularities

Alexandre Cornet; Alexandre Cornet

doi:10.3934/mbe.2018047

Mathematical Biosciences and Engineering

2018, Volume 15, Issue 5: 1055-1076. doi: 10.3934/mbe.2018047

Previous Article Next Article

Mathematical modelling of cardiac pulse wave reflections due to arterial irregularities

Alexandre Cornet ^{1,2
,
,}

1.
Imperial College London, South Kensington Campus, London SW72AZ, United Kingdom
2.
Ecole normale supérieure Paris-Saclay, 61 Avenue du Président Wilson, Cachan 94230, France

Received: 24 December 2016 Accepted: 23 January 2018 Published: 01 October 2018
MSC : Primary: 15A24, 76D05; Secondary: 65F30

This research aims to model cardiac pulse wave reflections due to the presence of arterial irregularities such as bifurcations, stiff arteries, stenoses or aneurysms. When an arterial pressure wave encounters an irregularity, a backward reflected wave travels upstream in the artery and a forward wave is transmitted downstream. The same process occurs at each subsequent irregularity, leading to the generation of multiple waves. An iterative algorithm is developed and applied to pathological scenarios to predict the pressure waveform of the reflected wave due to the presence of successive arterial irregularities. For an isolated stenosis, analysing the reflected pressure waveform gives information on its severity. The presence of a bifurcation after a stenosis tends do diminish the amplitude of the reflected wave, as bifurcations' reflection coefficients are relatively small compared to the ones of stenoses or aneurysms. In the case of two stenoses in series, local extrema are observed in the reflected pressure waveform which appears to be a characteristic of stenoses in series along an individual artery. Finally, we model a progressive change in stiffness in the vessel's wall and observe that the less the gradient stiffness is important, the weaker is the reflected wave.

Keywords:

Citation: Alexandre Cornet. Mathematical modelling of cardiac pulse wave reflections due to arterial irregularities[J]. Mathematical Biosciences and Engineering, 2018, 15(5): 1055-1076. doi: 10.3934/mbe.2018047

Related Papers:

[1]	Lela Dorel . Glucose level regulation via integral high-order sliding modes. Mathematical Biosciences and Engineering, 2011, 8(2): 549-560. doi: 10.3934/mbe.2011.8.549
[2]	Danilo T. Pérez-Rivera, Verónica L. Torres-Torres, Abraham E. Torres-Colón, Mayteé Cruz-Aponte . Development of a computational model of glucose toxicity in the progression of diabetes mellitus. Mathematical Biosciences and Engineering, 2016, 13(5): 1043-1058. doi: 10.3934/mbe.2016029
[3]	Lingmin Lin, Kailai Liu, Huan Feng, Jing Li, Hengle Chen, Tao Zhang, Boyun Xue, Jiarui Si . Glucose trajectory prediction by deep learning for personal home care of type 2 diabetes mellitus: modelling and applying. Mathematical Biosciences and Engineering, 2022, 19(10): 10096-10107. doi: 10.3934/mbe.2022472
[4]	Xingjia Li, Jinan Gu, Zedong Huang, Chen Ji, Shixi Tang . Hierarchical multiloop MPC scheme for robot manipulators with nonlinear disturbance observer. Mathematical Biosciences and Engineering, 2022, 19(12): 12601-12616. doi: 10.3934/mbe.2022588
[5]	Xiaohan Yang, Yinghao Cui, Zhanhang Yuan, Jie Hang . RISE-based adaptive control of electro-hydraulic servo system with uncertain compensation. Mathematical Biosciences and Engineering, 2023, 20(5): 9288-9304. doi: 10.3934/mbe.2023407
[6]	Kimberly Fessel, Jeffrey B. Gaither, Julie K. Bower, Trudy Gaillard, Kwame Osei, Grzegorz A. Rempała . Mathematical analysis of a model for glucose regulation. Mathematical Biosciences and Engineering, 2016, 13(1): 83-99. doi: 10.3934/mbe.2016.13.83
[7]	Balázs Csutak, Gábor Szederkényi . Robust control and data reconstruction for nonlinear epidemiological models using feedback linearization and state estimation. Mathematical Biosciences and Engineering, 2025, 22(1): 109-137. doi: 10.3934/mbe.2025006
[8]	Hugo Flores-Arguedas, Marcos A. Capistrán . Bayesian analysis of Glucose dynamics during the Oral Glucose Tolerance Test (OGTT). Mathematical Biosciences and Engineering, 2021, 18(4): 4628-4647. doi: 10.3934/mbe.2021235
[9]	Huazhong Yang, Wang Li, Maojin Tian, Yangfeng Ren . A personalized multitasking framework for real-time prediction of blood glucose levels in type 1 diabetes patients. Mathematical Biosciences and Engineering, 2024, 21(2): 2515-2541. doi: 10.3934/mbe.2024111
[10]	Xingjia Li, Jinan Gu, Zedong Huang, Wenbo Wang, Jing Li . Optimal design of model predictive controller based on transient search optimization applied to robotic manipulators. Mathematical Biosciences and Engineering, 2022, 19(9): 9371-9387. doi: 10.3934/mbe.2022436

Abstract

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 . 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.

Figure 1. Feedback and feedforward control flow diagram. For a detailed description, see the text.

DownLoad: Full-Size Img PowerPoint

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]:

$\begin{eqnarray} \frac{dg}{dt} & = & -m_1 g- h_a g+J , \end{eqnarray}$

(2.1)

$\begin{eqnarray} \frac{dh_a}{dt} & = & -m_2 h_a+m_3h , \end{eqnarray}$

(2.2)

$\begin{eqnarray} \frac{dh}{dt} & = & -m_4 h+u . \end{eqnarray}$

(2.3)

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:

Table 1. Values of parameters of the model (3.32)-(3.36).

Parameter	Value	Description
$m_1$	0.000296 (/min)	Glucose utilization rate constant
$m_2$	0.045 (/min)	Insulin effect rate constant
$m_3$	0.000012 (/(pmol/l)/min $^2$ )	Insulin effect rate constant
$m_4$	0.268 (/min)	Insulin metabolism constant
$G_b$	0.3218 (mM/min)	Basal glucose input
$\bar{g}$	6 (mM)	Glucose equilibrium
$\bar{h}$	200 (pmol/l)	Insulin equilibrium
$\bar{h}_a$	0.0533 (/min)	Insulin effect equilibrium
$\bar{u}$	53.63 (pmol/l/min)	Insulin input equilibrium

| Show Table

DownLoad: CSV

$\begin{equation} y = g . \end{equation}$

(2.4)

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:

$\begin{equation} J = G_b+G_d, \end{equation}$

(2.5)

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:

$\begin{eqnarray} 0 & = & -m_1 \bar{g}- \bar{g}\bar{h}_a +G_b , \end{eqnarray}$

(2.6)

$\begin{eqnarray} 0 & = & -m_2 \bar{h}_a+m_3\bar{h} , \end{eqnarray}$

(2.7)

$\begin{eqnarray} 0 & = & -m_4 \bar{h}+\bar{u} . \end{eqnarray}$

(2.8)

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):

$\begin{eqnarray} \frac{dg}{dt} & = & -(m_1+w_1) g- gh_a+ w_1 \bar{g} + G_b+G_d, \end{eqnarray}$

(2.9)

$\begin{eqnarray} \frac{dh_a}{dt} & = & - (m_2 +w_2) h_a+(m_3 +w_3) h+ w_2 \bar{h}_a-w_3 \bar{h} , \end{eqnarray}$

(2.10)

$\begin{eqnarray} \frac{dh}{dt} & = & - (m_4 +w_4) h+w_4 \bar{h}+u . \end{eqnarray}$

(2.11)

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$

$\lim\limits_{t\to \infty} (g(t)- r(t)) = 0.$

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

$\begin{eqnarray} G_d& = & a_0 +a_1\sin (\omega_d t) +b_1\cos(\omega_d t)+a_2\sin (2\omega_d t)+b_2\cos(2\omega_d t)+\cdots\\ &&+ a_n\sin (n\omega_d t)+b_n\cos(n\omega_d t). \end{eqnarray}$

(3.1)

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

$\frac{d\mathbf{v}_d}{dt} = \mathbf{A}_{d} \mathbf{v}_d,$

where

$\mathbf{A}_{d} = \left[\begin{array}{cccccccc} 0&0&0&0&0& \cdots& 0&0\\ 0& 0& \omega_d& 0 & 0& \cdots& 0&0 \\ 0&-\omega_d &0&0&0& \cdots& 0&0\\ 0&0&0&0&2\omega_d& \cdots& 0&0\\ 0&0&0&-2\omega_d&0& \cdots& 0&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&0&0& \cdots& 0&n\omega_d\\ 0&0&0&0&0& \cdots& -n\omega_d&0 \end{array} \right].$

Then $G_d$ can be expressed as

$G_d = \mathbf{C}_d\mathbf{v}_d,$

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:

$r = \bar{g}+\mathbf{C}_r\mathbf{v}_r,$

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

$\frac{d\mathbf{v}_r}{dt} = \mathbf{A}_{r} \mathbf{v}_r$

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

$\mathbf{A}_{r} = \left[\begin{array}{cccccccc} 0&0&0&0&0& \cdots& 0&0\\ 0& 0& \omega_r& 0 & 0& \cdots& 0&0 \\ 0&-\omega_r &0&0&0& \cdots& 0&0\\ 0&0&0&0&2\omega_r& \cdots& 0&0\\ 0&0&0&-2\omega_r&0& \cdots& 0&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&0&0& \cdots& 0&m\omega_r\\ 0&0&0&0&0& \cdots& -m\omega_r&0 \end{array} \right].$

Let

$\mathbf{v} = \left[\begin{array}{c} \mathbf{v}_r\\ \mathbf{v}_d \end{array} \right],\quad \mathbf{A}_e = \left[\begin{array}{cc} \mathbf{A}_{r}&0\\ 0& \mathbf{A}_{d} \end{array} \right].$

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

$\begin{equation} \frac{d\mathbf{v} }{dt} = \mathbf{A}_e \mathbf{v}. \end{equation}$

(3.2)

We introduce the tracking error:

$e = g-r = g-\bar{g}-\mathbf{C}_r\mathbf{v}_r.$

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

$\begin{eqnarray} \frac{dg}{dt} & = & -(m_1+w_1) g- gh_a+ w_1 \bar{g} + G_b+[\mathbf{0},\mathbf{C}_d]\mathbf{v}, \end{eqnarray}$

(3.3)

$\begin{eqnarray} \frac{dh_a}{dt} & = & - (m_2 +w_2) h_a+(m_3 +w_3) h+ w_2 \bar{h}_a-w_3 \bar{h} , \end{eqnarray}$

(3.4)

$\begin{eqnarray} \frac{dh}{dt} & = & - (m_4 +w_4) h+w_4 \bar{h}+u , \end{eqnarray}$

(3.5)

$\begin{eqnarray} \frac{d\mathbf{v} }{dt}& = &\mathbf{A}_{e} \mathbf{v}, \end{eqnarray}$

(3.6)

$\begin{eqnarray} e & = &g- \bar{g}-[\mathbf{C}_r, \mathbf{0}]\mathbf{v}. \end{eqnarray}$

(3.7)

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:

$\begin{equation} g = \hat{g}+\bar{g} +G(\mathbf{v}),\quad h_a = \hat{h}_a+\bar{h}_a+H_a(\mathbf{v}),\quad h = \hat{h}+\bar{h} +H(\mathbf{v}),\quad u = \hat{u}+\bar{u} +U( \mathbf{v}) \end{equation}$

(3.8)

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

$\begin{eqnarray} \frac{d\hat{g}}{dt}+\frac{\partial G}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v} & = & -a_1(\mathbf{w}) (\hat{g}+G(\mathbf{v}))- \bar{g} (\hat{h}_a+H_a(\mathbf{v}))\\ &&- [ \hat{g}\hat{h}+G(\mathbf{v}) \hat{h}_a+\hat{g}H_a(\mathbf{v})] - G(\mathbf{v})H_a(\mathbf{v})\\ && +[\mathbf{0},\mathbf{C}_d]\mathbf{v} , \end{eqnarray}$

(3.9)

$\begin{eqnarray} \frac{d\hat{h}_a}{dt}+\frac{\partial H_a}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v} & = & -a_2(\mathbf{w}) (\hat{h}_a+H_a(\mathbf{v}))+a_3(\mathbf{w}) (\hat{h}+H(\mathbf{v})) , \end{eqnarray}$

(3.10)

$\begin{eqnarray} \frac{d\hat{h}}{dt}+\frac{\partial H}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v} & = & -a_4(\mathbf{w}) (\hat{h}+H(\mathbf{v}))+\hat{u} +U(\mathbf{v}) , \end{eqnarray}$

(3.11)

$\begin{eqnarray} e & = &\hat{g}+G(\mathbf{v}) -[\mathbf{C}_r, \mathbf{0}]\mathbf{v} , \end{eqnarray}$

(3.12)

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

$\begin{eqnarray} a_1(\mathbf{w}) & = & m_1+w_1+ \bar{h}_a, \end{eqnarray}$

(3.13)

$\begin{eqnarray} a_2(\mathbf{w}) & = & m_2 +w_2 , \end{eqnarray}$

(3.14)

$\begin{eqnarray} a_3(\mathbf{w}) & = & m_3 +w_3 , \end{eqnarray}$

(3.15)

$\begin{eqnarray} a_4(\mathbf{w}) & = & m_4 +w_4. \end{eqnarray}$

(3.16)

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

$\begin{eqnarray} \frac{\partial G}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v}& = & -a_1(\mathbf{w}) G -\bar{g} H_a- GH_a + [\mathbf{0},\mathbf{C}_d]\mathbf{v} , \end{eqnarray}$

(3.17)

$\begin{eqnarray} \frac{\partial H_a}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v} & = & -a_2(\mathbf{w}) H_a+a_3(\mathbf{w}) H , \end{eqnarray}$

(3.18)

$\begin{eqnarray} \frac{\partial H}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v} & = & -a_4(\mathbf{w}) H+U , \end{eqnarray}$

(3.19)

$\begin{eqnarray} 0 & = &G-[\mathbf{C}_r, \mathbf{0}]\mathbf{v}. \end{eqnarray}$

(3.20)

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:

$\begin{eqnarray} \frac{d\hat{g}}{dt} & = & -a_1(\mathbf{w}) \hat{g} -\bar{g} \hat{h}_a- ( \hat{g}\hat{h}+G(\mathbf{v}) \hat{h}_a+\hat{g}H_a(\mathbf{v})), \end{eqnarray}$

(3.21)

$\begin{eqnarray} \frac{d\hat{h}_a}{dt} & = & -a_2(\mathbf{w}) \hat{h}_a + a_3(\mathbf{w}) \hat{h} , \end{eqnarray}$

(3.22)

$\begin{eqnarray} \frac{d\hat{h}}{dt} & = & -a_4(\mathbf{w}) \hat{h} + \hat{u} , \end{eqnarray}$

(3.23)

$\begin{eqnarray} e & = &\hat{g} . \end{eqnarray}$

(3.24)

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

$\begin{eqnarray} G(\mathbf{v}) & = & [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}, \end{eqnarray}$

(3.25)

$\begin{eqnarray} F(\mathbf{v}) & = & [\mathbf{0},\;\mathbf{C}_d]\mathbf{v}-a_1(\mathbf{w}) [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v} -[\mathbf{C}_r, \; \mathbf{0}] \mathbf{A}_e \mathbf{v} , \end{eqnarray}$

(3.26)

$\begin{eqnarray} H_a(\mathbf{v}) & = &\frac{ F(\mathbf{v}) } {\bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}} , \end{eqnarray}$

(3.27)

$\begin{eqnarray} H(\mathbf{v}) & = & \frac{ \frac{\partial H_a}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v}+a_2(\mathbf{w}) H_a }{a_3(\mathbf{w})} \end{eqnarray}$

(3.28)

$\begin{eqnarray} & = & \frac{ ( \bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}) \frac{\partial F}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v} -F(\mathbf{v})[\mathbf{C}_r, \; \mathbf{0}] \mathbf{A}_e \mathbf{v} } {a_3(\mathbf{w})[\bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}]^2} \\ &&+ \frac{ a_2(\mathbf{w})F(\mathbf{v})} {a_3(\mathbf{w})( \bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v})}, \end{eqnarray}$

(3.29)

$\begin{eqnarray} U(\mathbf{v}) & = & \frac{\partial H}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v}+a_4(\mathbf{w}) H \\ & = & \frac{ ( [\mathbf{C}_r, \; \mathbf{0}] \mathbf{A}_e \mathbf{v}) \frac{\partial F}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v} +( \bar{g}+ [\mathbf{C}_r, \; \mathbf{0}] \mathbf{v})\frac{\partial F}{\partial \mathbf{v}} \mathbf{A}_e^2 \mathbf{v} -\frac{\partial F}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v}[\mathbf{C}_r, \; \mathbf{0}] \mathbf{A}_e \mathbf{v}-F(\mathbf{v})[\mathbf{C}_r, \; \mathbf{0}] \mathbf{A}_e^2 \mathbf{v} } {a_3 (\mathbf{w})[\bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}]^2} \\ &&- \frac{ 2 [( \bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}) \frac{\partial F}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v} -F(\mathbf{v})[\mathbf{C}_r, \; \mathbf{0}] \mathbf{A}_e \mathbf{v} ][\bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}]a_3(\mathbf{w}) [\mathbf{C}_r, \; \mathbf{0}]\mathbf{A}_e\mathbf{v}} {a_3^2(\mathbf{w})[\bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}]^4} \\ &&+ \frac{ a_2(\mathbf{w})\frac{\partial F}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v} a_3 (\mathbf{w})( \bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}) - a_2(\mathbf{w})F(\mathbf{v})a_3 (\mathbf{w}) [\mathbf{C}_r, \; \mathbf{0}]\mathbf{A}_e\mathbf{v})} {a_3^2(\mathbf{w})( \bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v})^2}\\ &&+a_4(\mathbf{w})\frac{ ( \bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}) \frac{\partial F}{\partial \mathbf{v}} \mathbf{A}_e \mathbf{v} -F(\mathbf{v})[\mathbf{C}_r, \; \mathbf{0}] \mathbf{A}_e \mathbf{v} } {a_3(\mathbf{w})[\bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}]^2} \\ &&+ \frac{ a_2(\mathbf{w})a_4(\mathbf{w})F(\mathbf{v})} {a_3(\mathbf{w})( \bar{g}+ [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v})}. \end{eqnarray}$

(3.30)

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

$\hat{u} = k_g\hat{g}.$

Therefore, we have constructed a static output feedback

$\begin{eqnarray} u & = & \bar{u}+\hat{u}+U(\mathbf{v})\\ & = & \bar{u}+k_g\hat{g}+U(\mathbf{v})\\ & = & \bar{u}+k_g(g-\bar{g}- [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}))+U(\mathbf{v}) \end{eqnarray}$

(3.31)

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:

$\begin{eqnarray} \frac{dg}{dt} & = & -(m_1+w_1) g- gh_a+ w_1 \bar{g} + G_b+[\mathbf{0},\;\mathbf{C}_d]\mathbf{v}, \end{eqnarray}$

(3.32)

$\begin{eqnarray} \frac{dh_a}{dt} & = & - (m_2 +w_2) h_a+(m_3 +w_3) h+ w_2 \bar{h}_a-w_3 \bar{h} , \end{eqnarray}$

(3.33)

$\begin{eqnarray} \frac{dh}{dt} & = & - (m_4 +w_4) h+w_4 \bar{h}+ \bar{u}+k_g(g-\bar{g}- [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}))+U(\mathbf{v}) , \end{eqnarray}$

(3.34)

$\begin{eqnarray} \frac{d\mathbf{v} }{dt}& = &\mathbf{A}_{e} \mathbf{v}, \end{eqnarray}$

(3.35)

$\begin{eqnarray} e & = &g- \bar{g}- [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}. \end{eqnarray}$

(3.36)

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,

$\lim\limits_{t\to \infty} [g(t)-\bar{g}-[\mathbf{C}_r, \; \mathbf{0}] \mathbf{v}(t)] = 0.$

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

Figure 2. Glucose infusion shown in the second panel of Figure 2 of Sturis et al.'s ^[39]. The unit cc per hour is converted into mM per minute.

DownLoad: Full-Size Img PowerPoint

$J = 125[1+0.33\cdot\sin(\pi t/72 )]\; (\mbox{cc/hr}) .$

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

$J = \frac{125\cdot 0.2\cdot 1000}{60\cdot 18\cdot 60} [1+0.33\cdot\sin(\pi t/72 )] = 0.3858 [1+0.33\cdot\sin(\pi t/72 )]\; (\mbox{mM/min}) .$

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

$\mathbf{C}_d = [0.3858 - G_b, 0.1273, 0 ],$

and

$\mathbf{A}_{d} = \left[\begin{array}{ccc } 0&0&0 \\ 0& 0& \pi /72 \\ 0&-\pi /72 &0 \end{array} \right] .$

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

$g = a_0+\sum\limits_{i = 1}^{17}[a_i\sin(i\pi t/750 )+b_i\cos(i\pi t/750)],$

we obtain

$\begin{eqnarray*} a_0 & = & 7.599 , \\ a_1 & = & -0.2394 , \quad b_1 = -1.024,\\ a_2 & = & 0.328 , \quad b_2 = -0.2217,\\ a_3 & = & -0.2314 , \quad b_3 = -0.1684 ,\\ a_4 & = & -0.1663 , \quad b_4 = -0.05297,\\ a_5 & = & -0.1167 , \quad b_5 = 0.1937,\\ a_6 & = & 0.07061 , \quad b_6 = -0.1443,\\ a_7 & = & 0.3512 , \quad b_7 = -0.08337,\\ a_8 & = & 0.02829 , \quad b_8 = -0.0705,\\ a_9 & = & 0.2508 , \quad b_9 = -0.259,\\ a_{10} & = & 0.6841 , \quad b_{10} = -0.8653,\\ a_{11} & = & -0.3087 , \quad b_{11} = 0.5024,\\ a_{12} & = & -0.2117 , \quad b_{12} = 0.2057,\\ a_{13} & = & -0.07391 , \quad b_{13} = 0.02245,\\ a_{14} & = & -0.01738 , \quad b_{14} = -0.1047,\\ a_{15} & = & -0.1031 , \quad b_{15} = -0.05237,\\ a_{16} & = & -0.001116 , \quad b_{16} = 0.004296,\\ a_{17} & = & -0.04128 , \quad b_{17} = 0.003007. \end{eqnarray*}$

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

$\mathbf{C}_r = [a_0 - \bar{g}, a_1, b_1 , a_2, b_2 , \cdots a_{17}, b_{17} ],$

and

$\mathbf{A}_{r} = \left[\begin{array}{cccccccc} 0&0&0&0&0& \cdots& 0&0\\ 0& 0& \pi/750& 0 & 0& \cdots& 0&0 \\ 0&-\pi/750 &0&0&0& \cdots& 0&0\\ 0&0&0&0&2\pi/750& \cdots& 0&0\\ 0&0&0&-2\pi/750&0& \cdots& 0&0\\ \vdots&\vdots&\vdots&\vdots&\vdots&\ddots&\vdots&\vdots\\ 0&0&0&0&0& \cdots& 0&17\pi/750\\ 0&0&0&0&0& \cdots& -17\pi/750&0 \end{array} \right] .$

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

Figure 3. Fourier polynomial fitting of the blood glucose data shown in the second panel of Figure 2 of Sturis et al.'s ^[39]. The data were read by using the software Engauge Digitizer 4.1.

DownLoad: Full-Size Img PowerPoint

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 . 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. 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.

Figure 4. Glucose tracking under the static output feedback controller (3.31). The blood glucose

$g$ is tracking the desired time-varying reference determined by the experimental data from Sturis et al.'s ^[39]. The data were read by using the software Engauge Digitizer 4.1.

DownLoad: Full-Size Img PowerPoint

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), 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]).

Figure 5. Glucose fails to track its reference when the uncertainties are set to 0 in the static output feedback controller (3.31).

DownLoad: Full-Size Img PowerPoint

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

$\begin{eqnarray} \frac{d\tilde{g}}{dt} & = & -a_1(\mathbf{w}) \tilde{g} - \bar{g} \tilde{h}_a +[\mathbf{0},\;\mathbf{C}_d]\mathbf{v}, \end{eqnarray}$

(4.1)

$\begin{eqnarray} \frac{d\tilde{h}_a}{dt} & = & - a_2(\mathbf{w}) \tilde{h}_a + a_3(\mathbf{w}) \tilde{h} , \end{eqnarray}$

(4.2)

$\begin{eqnarray} \frac{d\tilde{h}}{dt} & = & - a_4(\mathbf{w}) \tilde{h} +\tilde{u} , \end{eqnarray}$

(4.3)

$\begin{eqnarray} \frac{d\mathbf{v} }{dt}& = &\mathbf{A}_{e} \mathbf{v}, \end{eqnarray}$

(4.4)

$\begin{eqnarray} e & = &\tilde{g}- [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}, \end{eqnarray}$

(4.5)

where

$\begin{equation} \tilde{g} = g-\bar{g},\; \tilde{h}_a = h_a-\bar{h}_a,\; \tilde{h} = h-\bar{h},\; \tilde{u} = u-\bar{u}. \end{equation}$

(4.6)

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

$\begin{equation} \frac{d\mathbf{z} }{dt} = \mathbf{G}_1 \mathbf{z}+\mathbf{G}_2(\tilde{g} - [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}), \end{equation}$

(4.7)

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:

$\begin{equation} \frac{d \mathbf{o}}{dt} = \mathbf{A} \mathbf{o}+ \mathbf{B}\tilde{u}+\mathbf{L}(\tilde{g}-\mathbf{C}\mathbf{o}) , \end{equation}$

(4.8)

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

$\begin{equation} \mathbf{A} = \left[\begin{array}{ccc} -a_1(\mathbf{0})& -\bar{g} &0\\ 0& -a_2(\mathbf{0}) &a_3(\mathbf{0})\\ 0& 0&-a_4(\mathbf{0}) \end{array} \right],\quad \mathbf{B} = \left[\begin{array}{c} 0\\ 0\\ 1 \end{array} \right],\quad \mathbf{C} = [1, 0, 0]. \end{equation}$

(4.9)

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

$\begin{equation} \tilde{u} = -\mathbf{K}_z \mathbf{z} -\mathbf{K}_o \mathbf{o}, \\ \end{equation}$

(4.10)

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

$\begin{eqnarray*} p(\lambda) & = &\lambda(\lambda^2+\omega_d^2)(\lambda^2+(2\omega_d)^2)\cdots (\lambda^2+(l\omega_d)^2)\\ & = &\lambda^{2l+1}+c_1\lambda^{2l}+\cdots+c_{2l}\lambda +c_{2l+1}, \end{eqnarray*}$

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

$\begin{eqnarray*} p(\lambda) & = &\lambda(\lambda^2+\omega_d^2)(\lambda^2+(2\omega_d)^2)\cdots (\lambda^2+(n\omega_d)^2)(\lambda^2+\omega_r^2)(\lambda^2+(2\omega_r)^2)\cdots (\lambda^2+(m\omega_r)^2)\\ & = &\lambda^{2l+1}+c_1\lambda^{2l}+\cdots+c_{2l}\lambda +c_{2l+1}, \end{eqnarray*}$

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:

$\begin{equation} \mathbf{G}_1 = \left[\begin{array}{ccccc} 0&1&0&\cdots&0\\ 0& 0& 1& \cdots & 0\\ \vdots&\vdots &\ddots&\vdots&\vdots\\ 0&0&0&\cdots&1\\ -c_{2l+1}&-c_{2l} &-c_{2l-1}&\cdots&-c_1 \end{array} \right],\quad \mathbf{G}_2 = \left[\begin{array}{c} 0\\ 0 \\ \vdots\\ 0\\ 1 \end{array} \right]. \end{equation}$

(4.11)

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

$\begin{equation} \tilde{g} = \hat{g} +\mathbf{G} \mathbf{v},\quad \tilde{h}_a = \hat{h}_a +\mathbf{H}_a \mathbf{v},\quad \tilde{h} = \hat{h} +\mathbf{H} \mathbf{v},\quad \quad \mathbf{z} = \hat{\mathbf{z}} +\mathbf{Z} \mathbf{v} ,\quad \quad \mathbf{o} = \hat{\mathbf{o}} +\mathbf{O} \mathbf{v}. \end{equation}$

(4.12)

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

$\begin{eqnarray*} \frac{d\hat{g}}{dt} + \mathbf{G} \mathbf{A}_{e}\mathbf{v}& = & - a_1(\mathbf{w}) (\hat{g} +\mathbf{G} \mathbf{v})- \bar{g} (\hat{h}_a +\mathbf{H}_a \mathbf{v}) +[\mathbf{0},\;\mathbf{C}_d]\mathbf{v}, \\ \frac{d\hat{h}_a}{dt}+ \mathbf{H}_a \mathbf{A}_{e}\mathbf{v} & = & - a_2(\mathbf{w}) (\hat{h}_a +\mathbf{H}_a \mathbf{v}) + a_3(\mathbf{w}) (\hat{h} +\mathbf{H} \mathbf{v}),\\ \frac{d\hat{h}}{dt}+ \mathbf{H} \mathbf{A}_{e}\mathbf{v} & = & - a_4(\mathbf{w}) (\hat{h} +\mathbf{H} \mathbf{v}) - \mathbf{K}_z (\hat{\mathbf{z}}+\mathbf{Z}\mathbf{v})- \mathbf{K}_o (\hat{\mathbf{o}}+\mathbf{O} \mathbf{v}) , \\ \frac{d\hat{\mathbf{z}}}{dt}+ \mathbf{Z} \mathbf{A}_{e}\mathbf{v} & = & \mathbf{G}_1 (\hat{\mathbf{z}}+\mathbf{Z} \mathbf{v})+ \mathbf{G}_2(\hat{g}+\mathbf{G} \mathbf{v} -[\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}),\\ \frac{d\hat{\mathbf{o}}}{dt}+ \mathbf{O} \mathbf{A}_{e}\mathbf{v} & = & \mathbf{A} (\hat{\mathbf{o}}+\mathbf{O} \mathbf{v})- \mathbf{B}[\mathbf{K}_z (\hat{\mathbf{z}}+\mathbf{Z}\mathbf{v})+ \mathbf{K}_o (\hat{\mathbf{o}}+\mathbf{O} \mathbf{v}) ]+\mathbf{L}(\hat{g} +\mathbf{G} \mathbf{v}-\mathbf{C}(\hat{\mathbf{o}}+\mathbf{O} \mathbf{v})),\\ e & = &\hat{g} +\mathbf{G} \mathbf{v} - [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}. \end{eqnarray*}$

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

$\begin{eqnarray} \mathbf{G} \mathbf{A}_{e} & = & - a_1(\mathbf{w}) \mathbf{G} - \bar{g} \mathbf{H}_a +[\mathbf{0},\;\mathbf{C}_d], \end{eqnarray}$

(4.13)

$\begin{eqnarray} \mathbf{H}_a \mathbf{A}_{e} & = & - a_2(\mathbf{w}) \mathbf{H}_a +a_3(\mathbf{w}) \mathbf{H} , \end{eqnarray}$

(4.14)

$\begin{eqnarray} \mathbf{H} \mathbf{A}_{e} & = & - a_4(\mathbf{w}) \mathbf{H} -\mathbf{K}_z \mathbf{Z} -\mathbf{K}_o \mathbf{O}, \end{eqnarray}$

(4.15)

$\begin{eqnarray} \mathbf{Z} \mathbf{A}_{e} & = & \mathbf{G}_1 \mathbf{Z} + \mathbf{G}_2( \mathbf{G} -[\mathbf{C}_r, \; \mathbf{0}] ), \end{eqnarray}$

(4.16)

$\begin{eqnarray} \mathbf{O} \mathbf{A}_{e} & = & \mathbf{A} \mathbf{O} -\mathbf{B}\mathbf{K}_z \mathbf{Z} -\mathbf{B}\mathbf{K}_o \mathbf{O}+\mathbf{L}(\mathbf{G}-\mathbf{C}\mathbf{O}), \end{eqnarray}$

(4.17)

$\begin{eqnarray} 0 & = & \mathbf{G} - [\mathbf{C}_r, \; \mathbf{0}] \end{eqnarray}$

(4.18)

and then

$\begin{eqnarray} \frac{d\hat{g}}{dt} & = & - a_1(\mathbf{w}) \hat{g}- \bar{g} \hat{h}_a , \end{eqnarray}$

(4.19)

$\begin{eqnarray} \frac{d\hat{h}_a}{dt} & = & - a_2(\mathbf{w}) \hat{h}_a + a_3(\mathbf{w}) \hat{h} , \end{eqnarray}$

(4.20)

$\begin{eqnarray} \frac{d\hat{h}}{dt} & = & - a_4(\mathbf{w}) \hat{h} - \mathbf{K}_z \hat{\mathbf{z}} - \mathbf{K}_o \hat{\mathbf{o}} , \end{eqnarray}$

(4.21)

$\begin{eqnarray} \frac{d\hat{\mathbf{z}}}{dt} & = & \mathbf{G}_1 \hat{\mathbf{z}} + \mathbf{G}_2 \hat{g} , \end{eqnarray}$

(4.22)

$\begin{eqnarray} \frac{d\hat{\mathbf{o}}}{dt} & = & \mathbf{A} \hat{\mathbf{o}}- \mathbf{B}[\mathbf{K}_z \hat{\mathbf{z}}+ \mathbf{K}_o \hat{\mathbf{o}} ]+\mathbf{L}(\hat{g} -\mathbf{C}\hat{\mathbf{o}}), \end{eqnarray}$

(4.23)

$\begin{eqnarray} e & = &\hat{g} . \end{eqnarray}$

(4.24)

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:

$\begin{eqnarray} \frac{d\hat{\mathbf{x}}}{dt} & = & \mathbf{A} \hat{\mathbf{x}}- \mathbf{B}[\mathbf{K}_z \hat{\mathbf{z}}+ \mathbf{K}_o \hat{\mathbf{o}} ], \end{eqnarray}$

(4.25)

$\begin{eqnarray} \frac{d\hat{\mathbf{z}}}{dt} & = & \mathbf{G}_1 \hat{\mathbf{z}} + \mathbf{G}_2 \mathbf{C}\hat{\mathbf{x}} , \end{eqnarray}$

(4.26)

$\begin{eqnarray} \frac{d\hat{\mathbf{o}}}{dt} & = & \mathbf{A} \hat{\mathbf{o}}- \mathbf{B}[\mathbf{K}_z \hat{\mathbf{z}}+ \mathbf{K}_o \hat{\mathbf{o}} ]+\mathbf{L}(\mathbf{C}\hat{\mathbf{x}} -\mathbf{C}\hat{\mathbf{o}}). \end{eqnarray}$

(1.27)

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

$\begin{eqnarray} \frac{d\hat{\mathbf{x}}}{dt} & = & ( \mathbf{A} - \mathbf{B}\mathbf{K}_o)\hat{\mathbf{x}}- \mathbf{B}[\mathbf{K}_z \hat{\mathbf{z}}- \mathbf{K}_o \hat{\mathbf{q}} ], \end{eqnarray}$

(4.28)

$\begin{eqnarray} \frac{d\hat{\mathbf{z}}}{dt} & = & \mathbf{G}_1 \hat{\mathbf{z}} + \mathbf{G}_2 \mathbf{C}\hat{\mathbf{x}} , \end{eqnarray}$

(4.29)

$\begin{eqnarray} \frac{d\hat{\mathbf{q}}}{dt} & = & (\mathbf{A} -\mathbf{L}\mathbf{C})\hat{\mathbf{q}}, \end{eqnarray}$

(4.30)

and then

$\begin{equation} \frac{d}{dt}\left[\begin{array}{c} \hat{\mathbf{x}} \\ \hat{\mathbf{z}}\\ \hat{\mathbf{q}} \end{array} \right] = \left[\begin{array}{ccc} \mathbf{A} - \mathbf{B}\mathbf{K}_o& - \mathbf{B}\mathbf{K}_z & \mathbf{B}\mathbf{K}_o\\ \mathbf{G}_2 \mathbf{C}&\mathbf{G}_1 &\mathbf{0}\\ \mathbf{0}&\mathbf{0} &\mathbf{A} -\mathbf{L}\mathbf{C} \end{array} \right]\left[\begin{array}{c} \hat{\mathbf{x}} \\ \hat{\mathbf{z}}\\ \hat{\mathbf{q}} \end{array} \right]. \end{equation}$

(4.31)

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

$\begin{equation} \mathbf{\Phi} = \left[\begin{array}{cc} \mathbf{A} & \mathbf{0} \\ \mathbf{G}_2 \mathbf{C}&\mathbf{G}_1 \end{array} \right],\quad \mathbf{\Psi} = \left[\begin{array}{ccc} \mathbf{B} \\ \mathbf{0} \end{array} \right] \end{equation}$

(4.32)

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

$\begin{equation} \mathbf{\Omega} = \left[\begin{array}{ccc} \mathbf{A} - \mathbf{B}\mathbf{K}_o& - \mathbf{B}\mathbf{K}_z & \mathbf{B}\mathbf{K}_o\\ \mathbf{G}_2 \mathbf{C}&\mathbf{G}_1 &\mathbf{0}\\ \mathbf{0}&\mathbf{0} &\mathbf{A}-\mathbf{L}\mathbf{C} \end{array} \right] \end{equation}$

(4.33)

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:

$\begin{equation} \left[\begin{array}{c } \mathbf{X} \\ \mathbf{Z}\\ \mathbf{Z} \end{array} \right] \mathbf{A}_e = \mathbf{Q}\left[\begin{array}{c } \mathbf{X} \\ \mathbf{Z} \\ \mathbf{O} \end{array} \right]+ \left[\begin{array}{c } \mathbf{C}^T [0,\;\mathbf{C}_d] \\ \mathbf{G}_2[\mathbf{C}_r, \; 0]\\ \mathbf{0} \end{array} \right], \end{equation}$

(4.34)

where

$\begin{equation} \mathbf{Q} = \left[\begin{array}{ccc} \mathbf{A} & - \mathbf{B}\mathbf{K}_z & -\mathbf{B}\mathbf{K}_o\\ \mathbf{G}_2 \mathbf{C}&\mathbf{G}_1 &\mathbf{0}\\ \mathbf{L}\mathbf{C} &- \mathbf{B}\mathbf{K}_z &\mathbf{A} -\mathbf{L}\mathbf{C}-\mathbf{B}\mathbf{K}_o \end{array} \right]. \end{equation}$

(4.35)

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

$\begin{eqnarray} \mathbf{Z}_{2} & = &\mathbf{Z}_{1} \mathbf{A}_e, \end{eqnarray}$

(4.36)

$\begin{eqnarray} \mathbf{Z}_{3} & = &\mathbf{Z}_{2} \mathbf{A}_e = \mathbf{Z}_{1} \mathbf{A}_e^2, \end{eqnarray}$

(4.37)

$\begin{eqnarray} &\vdots& \\ \mathbf{Z}_{2l+1} & = &\mathbf{Z}_{2l} \mathbf{A}_e = \mathbf{Z}_{1} \mathbf{A}_e^{2l}, \end{eqnarray}$

(4.38)

$\begin{eqnarray} \mathbf{G} & = &[\mathbf{C}_r, \; 0]+\mathbf{Z}_{2l+1} \mathbf{A}_e+c_{2l+1}\mathbf{Z}_{1}+c_{2l}\mathbf{Z}_{2}+\cdots+c_1\mathbf{Z}_{2l+1} \\ & = &[\mathbf{C}_r, \; 0]+\mathbf{Z}_{1} \mathbf{A}_e^{2l+1} +c_{2l+1}\mathbf{Z}_{1}+c_{2l}\mathbf{Z}_{1}\mathbf{A}_e+\cdots +c_1\mathbf{Z}_{1} \mathbf{A}_e^{2l} \\ & = &[\mathbf{C}_r, \; 0]+\mathbf{Z}_{1} ( c_{2l+1}\mathbf{I}+c_{2l} \mathbf{A}_e+\cdots +c_1 \mathbf{A}_e^{2l}+\mathbf{A}_e^{2l+1}) \\ & = &[\mathbf{C}_r, \; 0]+\mathbf{Z}_{1}\mathbf{0} \\ & = &[\mathbf{C}_r, \; 0] . \end{eqnarray}$

(4.39)

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):

$\begin{eqnarray} u& = &\bar{u}+\tilde{u}\\ & = & \bar{u}-\mathbf{K}_z \mathbf{z}-\mathbf{K}_o \mathbf{o}, \end{eqnarray}$

(4.40)

$\begin{eqnarray} \frac{d\mathbf{z} }{dt} & = &\mathbf{G}_1 \mathbf{z}+\mathbf{G}_2( g-\bar{g} - [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}), \end{eqnarray}$

(4.41)

$\begin{eqnarray} \frac{d \mathbf{o}}{dt} & = & \mathbf{A} \mathbf{o}- \mathbf{B}(\mathbf{K}_z \mathbf{z}+\mathbf{K}_o \mathbf{o})+\mathbf{L}(g-\bar{g}-\mathbf{C}\mathbf{o}). \end{eqnarray}$

(4.42)

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

$\begin{eqnarray} \frac{dg}{dt} & = & -(m_1+w_1) g- gh_a+ w_1 \bar{g} +G_b+[\mathbf{0},\;\mathbf{C}_d]\mathbf{v}, \end{eqnarray}$

(4.43)

$\begin{eqnarray} \frac{dh_a}{dt} & = & - (m_2 +w_2) h_a+ (m_3 +w_3) h +w_2 \bar{h}_a-w_3 \bar{h} , \end{eqnarray}$

(4.44)

$\begin{eqnarray} \frac{dh}{dt} & = & - (m_4 +w_4) h+w_4 \bar{h}+\bar{u}-\mathbf{K}_z \mathbf{z}-\mathbf{K}_o \mathbf{o}, \end{eqnarray}$

(4.45)

$\begin{eqnarray} \frac{d\mathbf{z} }{dt} & = &\mathbf{G}_1 \mathbf{z}+\mathbf{G}_2( g-\bar{g} - [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}), \end{eqnarray}$

(4.46)

$\begin{eqnarray} \frac{d \mathbf{o}}{dt} & = & \mathbf{A} \mathbf{o}- \mathbf{B}(\mathbf{K}_z \mathbf{z}+\mathbf{K}_o \mathbf{o})+\mathbf{L}(g-\bar{g}-\mathbf{C}\mathbf{o}), \end{eqnarray}$

(4.47)

$\begin{eqnarray} \frac{d\mathbf{v} }{dt}& = &\mathbf{A}_{e} \mathbf{v}, \end{eqnarray}$

(4.48)

$\begin{eqnarray} e & = &g- \bar{g}- [\mathbf{C}_r, \; \mathbf{0}]\mathbf{v}. \end{eqnarray}$

(4.49)

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

$\begin{equation*} \mathbf{\Omega} = \left[\begin{array}{ccc} \mathbf{A} - \mathbf{B}\mathbf{K}_o& - \mathbf{B}\mathbf{K}_z & \mathbf{B}\mathbf{K}_o\\ \mathbf{G}_2 \mathbf{C}&\mathbf{G}_1 &\mathbf{0}\\ \mathbf{0}&\mathbf{0} &\mathbf{A} -\mathbf{L}\mathbf{C} \end{array} \right] \end{equation*}$

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,

$\begin{equation*} \lim\limits_{t\to \infty} [g(t)-\bar{g}- [\mathbf{C}_r, \; \mathbf{0}] \mathbf{v}(t)] = 0. \end{equation*}$

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

$\begin{equation*} \mathbf{A}_{e} = \left[\begin{array}{ccccccccc} 0&0&0& \cdots& 0&0&0&0&0\\ 0& 0& \pi/750& \cdots& 0&0&0&0&0 \\ 0&-\pi/750 &0& \cdots& 0&0&0&0&0\\ 0&0&0& \cdots& 0&0&0&0&0\\ 0&0&0& \cdots& 0&0&0&0&0\\ \vdots&\vdots&\vdots& \ddots&\vdots&\vdots&\vdots&\vdots&\vdots\\ 0&0&0& \cdots& 0&17\pi/750&0&0&0\\ 0&0&0& \cdots& -17\pi/750&0&0&0&0\\ 0&0&0& \cdots& 0&0&0&0&0\\ 0&0&0& \cdots& 0&0&0&0&\pi/72\\ 0&0&0& \cdots& 0&0&0&-\pi/72&0 \end{array} \right] . \end{equation*}$

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

$\begin{equation*} \tilde{\mathbf{A}}_{e} = \left[\begin{array}{cccccccccc} 0& 0& 0& 0& 0& 0&0&0&0&0\\ 0& 0& \pi/750& 0& 0& 0&0&0&0&0 \\ 0&-\pi/750 &0& 0& 0& 0&0&0&0&0\\ 0& 0& 0& 0& 2\pi/750&0&0&0&0&0\\ 0& 0& 0& - 2\pi/750&0 & 0&0&0&0&0\\ 0& 0& 0& 0& 0&0&3\pi/750&0&0&0\\ 0& 0& 0& 0& 0&-3\pi/750&0&0&0&0\\ 0&0&0& 0& 0&0&0&0&0&0\\ 0&0&0& 0& 0&0&0&0&0&\pi/72\\ 0&0&0& 0& 0&0&0&0&-\pi/72&0 \end{array} \right] . \end{equation*}$

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

$\begin{eqnarray*} \mathbf{G}_1 & = & \left[\begin{array}{ccccccccc} 0& 1& 0& 0& 0& 0&0&0&0\\ 0& 0& 1& 0& 0& 0&0&0&0\\ 0& 0& 0& 1& 0& 0&0&0&0\\ 0& 0& 0& 0& 1& 0&0&0&0\\ 0& 0& 0& 0& 0& 1&0&0&0\\ 0& 0& 0& 0& 0& 0&1 &0&0\\ 0& 0& 0& 0& 0& 0&0&1&0\\ 0& 0& 0& 0& 0& 0&0&0&1\\ 0 &-3.70\cdot 10^{-16} & 0& -2.89\cdot 10^{-11} & 0 &-4.83\cdot 10^{-7} & 0 & -0.0021495 & 0 \end{array} \right] ,\\ \mathbf{G}_2 & = & \left[\begin{array}{c} 0\\ 0 \\ 0\\ 0\\ 0\\ 0 \\ 0\\ 0\\ 1 \end{array} \right]. \end{eqnarray*}$

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

$\begin{eqnarray*} \mathbf{K}_o & = & [ -55379.84, 379273.25, 2.93],\\ \mathbf{K}_z & = & [ -0.00028, -0.019, -0.56, -9.41, -100.48, -728.94, -3701.22, -13180, -33587.83]. \end{eqnarray*}$

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

$\begin{equation*} \mathbf{L } = [ -0.30679, -0.011386,211.93]^T. \end{equation*}$

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).

Figure 6. Glucose tracking under the robust controller (4.40). The dynamics of the blood glucose tracking is almost identical to the one under the static controller (3.31).

DownLoad: Full-Size Img PowerPoint

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.

References

[1]	[ L. Augsburger,P. Reymon,E. Fonck,Z. Kulcar,M. Farhat,M. Ohta,N. Stergiopulos,D. A. Rufenacht, Methodologies to assess blood flow in cerebral aneurysms: Current state of research and perspectives, Journal of Neuroradiology, 36 (2009): 270-277.
[2]	[ I. Bakirtas,A. Antar, Effect of stenosis on solitary waves in arteries, International Journal of Engineering Science, 43 (2005): 730-743.
[3]	[ W.-S. Duan,Y.-R. Shi,X.-R Hong,K.-P. Lu,J.-B. Zhao, The reflection of soliton at multi-arterial bifurcations and the effect of the arterial inhomogeneity, Physics Letters A, 295 (2002): 133-138.
[4]	[ L. Formaggia,F. Nobile,A. Quarteroni,A. Veneziani, Multiscale modelling of the circulatory system: A preliminary analysis, Comput. Visual. Sci., 2 (1999): 75-83.
[5]	[ L. Formaggia,F. Nobile,A. Quarteroni,J.-F. Gerbeau, On the coupling of 3D and 1D Navier-Stokes equations for flow problems in compliant vessels, Comp. Methods Appl. Mech. Engng., 191 (2001): 561-582.
[6]	[ K. Hayashi, K. Handa, S. Nagasawa and A. Okumura, Stiffness and elastic behaviour of human intracranial and extracranial arteries, J. Biomech., 13 (1980), 175-179,181-184.
[7]	[ G. L. Langewouters,K. H. Wesseling,W. J. A. Goedhard, The static properties of 45 human thoracic and 20 abdominal aortas in vitro and the parameters of a new model, J. Biomech., 17 (1984): 425-435.
[8]	[ C. A. D. Leguy,E. M. H. Bosboom,H. Gelderblom,A. P. G. Hoeks,F. N. Van de Vosse, Estimation of distributed arterial mechanical properties using a wave propagation model in a reverse way, Medical Engineering & Physics, 32 (2010): 957-967.
[9]	[ K. S. Matthus,J. Alastruey,J. Peiro,A. W. Khir,P. Segers,R. P. Verdonck,K. H. Parker,S. J. Sherwin, Pulse wave propagation in a model human arterial network: Assessment of 1-D visco-elastic simulations against in vitro measurements, J. Biomechanics, 44 (2011): 2250-2258.
[10]	[ H. G. Morales,I. Larrabide,A. J. Geers,M. L. Aguilar,A. F. Frangi, Newtonian and non-Newtonian blood flow in coiled cerebral aneurysms, J. Biomechanics, 46 (2013): 2158-2164.
[11]	[ W. W. Nichols, J. W. Petersen, S. J. Denardo and D. D. Christou Arterial stiffness, wave reflection amplitude and left ventricular afterload are increased in overweight individuals, Artery Research, 7 (2013), 222-229.
[12]	[ Z. Ovadia-Blechman,S. Einav,U. Zaretsky,D. Castel,E. Eldar, Characterization of arterial stenosis and elasticity by analysis of high-frequency pressure wave components, Computer in Biology and Medicine, 33 (2003): 375-393.
[13]	[ C. S. Park,S. J. Payne, Nonlinear and viscous effects on wave propagation in an elastic axisymmetric vessel, J. of Fluids and Structures, 27 (2011): 134-144.
[14]	[ K. H. Parker, An introduction to wave intensity analysis, Medical & Biological Engineering & Computing, 47 (2009): 175-199.
[15]	[ T. J. Pedley, Nonlinear pulse wave reflection at an arterial stenosis, J. of Biomechanical Engineering, 105 (1983): 353-359.
[16]	[ S. I. S. Pinto,E. Doutel,J. B. L. M. Campos,J. M. Miranda, Blood analog fluid flow in vessels with stenosis: Development of an Openfoam code to stimulate pulsatile flow and elasticity of the fluid, APCBEE Procedia, 7 (2013): 73-79.
[17]	[ A. Quarteroni,L. Formaggia, Mathematical modelling and numerical simulation of the cardiovascular system, Handbook of Numerical Analysis, 12 (2004): 3-127.
[18]	[ P. Segers,J. Kips,B. Trachet,A. Swillens,S. Vermeersch,D. Mahieu,E. Rietzschel,M. D. Buyzere,L. V. Bortel, Limitations and pitfalls of non-invasive measurement of arterial pressure wave reflections and pulse wave velocity, Artery Research, 3 (2009): 79-88.
[19]	[ D. Shahmirzadi,E. E. Konofagou, Quantification of arterial wall inhomogeneity size, distribution, and modulus contrast using FSI numerical pulse wave propagation, Artery Research, 8 (2014): 57-65.
[20]	[ N. Stergiopulos,D. F. Young,T. R. Rogge, Computer simulation of arterial flow with applications to arterial and aortic stenoses, J. Biomechanics, 25 (1992): 1477-1488.
[21]	[ N. Stergiopulos,F. Spiridon,F. Pythoud,J. J. Meister, On the wave transmission and reflection properties of stenoses, J. Biomechanics, 29 (1996): 31-38.
[22]	[ A. Swillens,P. Segers, Assessment of arterial pressure wave reflection: Methodological considerations, Artery Research, 2 (2008): 122-131.
[23]	[ A. Tozeren, Elastic properties of arteries and their influence on the cardiovascular system, J. Biomech. Eng., 106 (1984): 182-185.
[24]	[ C. Tu,M. Deville,L. Dheur,L. Vanderschuren, Finite element simulation of pulsatile flow through arterial stenosis, J. Biomechanics, 25 (1992): 1141-1152.
[25]	[ J. J. Wang,K. H. Parker, Wave propagation in a model of the arterial circulation, J. Biomechanics, 37 (2004): 457-470.

This article has been cited by:

1.	J. Sooknanan, D. M. G. Comissiong, Trending on Social Media: Integrating Social Media into Infectious Disease Dynamics, 2020, 82, 0092-8240, 10.1007/s11538-020-00757-4
2.	Aníbal Coronel, Luis Friz, Ian Hess, María Zegarra, On the existence and uniqueness of an inverse problem in epidemiology, 2021, 100, 0003-6811, 513, 10.1080/00036811.2019.1608964
3.	Winfried Just, Joan Saldaña, Ying Xin, Oscillations in epidemic models with spread of awareness, 2018, 76, 0303-6812, 1027, 10.1007/s00285-017-1166-x
4.	Fahad Al Basir, Santanu Ray, Ezio Venturino, Role of media coverage and delay in controlling infectious diseases: A mathematical model, 2018, 337, 00963003, 372, 10.1016/j.amc.2018.05.042
5.	Chao Liu, Jane Heffernan, Stochastic dynamics in a delayed epidemic system with Markovian switching and media coverage, 2020, 2020, 1687-1847, 10.1186/s13662-020-02894-5
6.	Fulgensia Kamugisha Mbabazi, J. Y. T. Mugisha, Mark Kimathi, Global Stability of Pneumococcal Pneumonia with Awareness and Saturated Treatment, 2020, 2020, 1110-757X, 1, 10.1155/2020/3243957
7.	M. A. Aziz-Alaoui, Jean M.-S. Lubuma, Berge Tsanou, Prevalence-based modeling approach of schistosomiasis: global stability analysis and integrated control assessment, 2021, 40, 2238-3603, 10.1007/s40314-021-01414-9
8.	David Juher, David Rojas, Joan Saldaña, Robustness of behaviorally induced oscillations in epidemic models under a low rate of imported cases, 2020, 102, 2470-0045, 10.1103/PhysRevE.102.052301
9.	Tanuja Das, Prashant K. Srivastava, Anuj Kumar, Nonlinear dynamical behavior of an SEIR mathematical model: Effect of information and saturated treatment, 2021, 31, 1054-1500, 043104, 10.1063/5.0039048
10.	Shikha Jain, Sachin Kumar, Dynamic analysis of the role of innate immunity in SEIS epidemic model, 2021, 136, 2190-5444, 10.1140/epjp/s13360-021-01390-3
11.	Ellen P. Carlin, Koya C. Allen, Jeffrey J. Morgan, Jean-Paul Chretien, Suzan Murray, Deborah Winslow, Dawn Zimmerman, Behavioral Risk Modeling for Pandemics: Overcoming Challenges and Advancing the Science, 2021, 19, 2326-5094, 447, 10.1089/hs.2020.0209
12.	Khadija Akdim, Adil Ez-Zetouni, Mehdi Zahid, A stochastic vaccinated epidemic model incorporating Lévy processes with a general awareness-induced incidence, 2021, 14, 1793-5245, 2150044, 10.1142/S1793524521500443
13.	Talal Daghriri, Michael Proctor, Sarah Matthews, Evolution of Select Epidemiological Modeling and the Rise of Population Sentiment Analysis: A Literature Review and COVID-19 Sentiment Illustration, 2022, 19, 1660-4601, 3230, 10.3390/ijerph19063230
14.	Ze Chen, Yuan Wang, Yanjun Guan, Michael Jie Guo, Rong Xu, Long‐term effect of childhood pandemic experience on medical major choice: Evidence from the 2003 severe acute respiratory syndrome outbreak in China, 2023, 1057-9230, 10.1002/hec.4659
15.	Tanuja Das, Prashant K. Srivastava, Effect of a novel generalized incidence rate function in SIR model: Stability switches and bifurcations, 2023, 166, 09600779, 112967, 10.1016/j.chaos.2022.112967
16.	Rabindra Kumar Gupta, Rajanish Kumar Rai, Pankaj Kumar Tiwari, Arvind Kumar Misra, Maia Martcheva, A mathematical model for the impact of disinfectants on the control of bacterial diseases, 2023, 17, 1751-3758, 10.1080/17513758.2023.2206859
17.	A. K. Misra, Kuldeep Kumar Tripathi, Maitri Verma, Modeling the effects of psychological fear and media-induced awareness on the dynamics of infectious diseases, 2023, 0, 1531-3492, 0, 10.3934/dcdsb.2023118
18.	Jin Yang, Likun Guan, Zhuo Chen, Yuanshun Tan, Zijian Liu, Robert A. Cheke, Bifurcation analysis of a nonlinear pulse SIR model with media coverage, 2023, 111, 0924-090X, 19543, 10.1007/s11071-023-08869-x
19.	Aníbal Coronel, Fernando Huancas, Ian Hess, Alex Tello, The diffusion identification in a SIS reaction-diffusion system, 2023, 21, 1551-0018, 562, 10.3934/mbe.2024024
20.	Makayla Preston, Alexandria Carter, Eric Numfor, Modeling the Effects of Media Awareness on SARS-CoV-2 Transmission in Georgia, 2024, 10, 2349-5103, 10.1007/s40819-024-01759-9
21.	Arjun Kumar, Uma S. Dubey, Balram Dubey, Modeling the effect of vaccinations, hospital beds, and treatments on the dynamics of infectious disease with delayed optimal control and sensitivity analysis, 2024, 139, 2190-5444, 10.1140/epjp/s13360-024-05508-1
22.	Saduri Das, Prashant K. Srivastava, Pankaj Biswas, The role of social media in a tuberculosis compartmental model: Exploring Hopf-bifurcation and nonlinear oscillations, 2025, 230, 03784754, 323, 10.1016/j.matcom.2024.11.015
23.	Pritam Saha, Kalyan Kumar Pal, Uttam Ghosh, Pankaj Kumar Tiwari, Dynamic analysis of deterministic and stochastic SEIR models incorporating the Ornstein–Uhlenbeck process, 2025, 35, 1054-1500, 10.1063/5.0243656
24.	MAITRI VERMA, KULDEEP KUMAR TRIPATHI, ARVIND KUMAR MISRA, MODELING THE EFFECTS OF PSYCHOLOGICAL FEAR, TV AND SOCIAL MEDIA ADVERTISEMENTS ON THE CONTROL OF INFECTIOUS DISEASES, 2025, 33, 0218-3390, 553, 10.1142/S0218339025500147
25.	Xinghua Chang, Jianrong Wang, Maoxing Liu, Xue Yan, An SIS infectious disease model with nonlinear incidence and disease awareness on complex networks, 2025, 196, 09600779, 116349, 10.1016/j.chaos.2025.116349

Reader Comments

Your name:*

Email:*
© 2018 the Author(s), licensee AIMS Press. This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0)