ABSTRACT

Mathematical modelling of glucose-insulin system is very important in medicine as a necessary tool to understand the homeostatic control of human body. It can also be used to design clinical trials and in the evaluation of the diabetes prevention. In the last three decades so much work has been done in this direction. One of the most notable models is the global six compartment-mathematical model with 22 ordinary differential equations due to John Thomas Sorensen. This paper proposes a more simplified three compartment-mathematical model with only 6 ordinary differential equations by introducing a tissue compartment comprising kidney, gut, brain and periphery. For model parameter identification, we use inverse problems technique to solve a specific optimal control problem where data are obtained by solving the global model of John Thomas Sorensen. Numerical results show that the proposed model is adaptable to data and can be used to adjust diabetes mellitus type I or type II for diabetic patients.

Keywords:

Glucose, Insulin, Mathematical Model, Parameter, Identification

1. Introduction

It is common knowledge that lifestyle factors largely influence our health. These factors are diet, physical activity, smoking and psychological stress. The lifestyle changes have the influence on metabolism of some systems of human body including glucose-insulin system, for example disordered glucoregulation. Therefore two variables that have a bearing on glucose homeostasis are affected, those are: pancreas beta cell response to glucose and sensitivity of body to insulin. The key organs that control blood glucose are pancreas and liver. The key hormones are insulin and glucagon. Large-scale clinical trials have demonstrated the benefits of tight control of glucose-insulin system, minimizing disease complications and improving quality of life [1] . Human body needs to maintain glucose concentration level in a narrow range. According to the World Health Organization or the International Diabetes Federation the cut-offs for fasting blood sugars, hypoglycemia is below 70 mg/dl and hyperglycemia is more than 126 mg/dl. The upper normal limit is 200mg/dl. Many researches are motivated by the large population of diabetes patients in the world and the big health expenses to study the glucose-insulin endocrine metabolic regulatory system [2] [3] [4] [5] [6] . They are interested in what cause the dysfunctions of the system [7] , how to detect the onset of either type of diabetes including the so called prediabetes [8] [9] [10] [11] , and eventually provide more reasonable, more effective, more efficient and cost-effective treatments to diabetes. For example, diabetes mellitus is a metabolic disorder caused by insufficient insulin production in islet cells in the pancreas or by tissue resistance against secreted insulin, which leads to excessive glucose concentration in the blood.

Since the 1960s, mathematical models have been developed to describe glucose-insulin dynamics [7] . Starting from the model proposed by Bolie in 1961 [12] mathematical modeling of glucose-insulin interaction in normal body has been studied. Some of these mathematical models are interested in analysing the glucose disappearance and insulin sensitivity during an intravenous glucose tolerance test [13] and capture of plasma glucose and insulin dynamics during, as well as after, periods of mild-to-moderate exercise [14] . Others focus on overestimation of glucose effectiveness and the underestimation of insulin sensitivity [15] , capture absorption, distribution and disposal dynamics [16] . More detailed compartmental models have been proposed by Cobelli et al. [17] , Hovorka et al. [18] and Sorensen [19] . Most of these models have been used for diabetes mellitus modeling by adjusting the model for type I or type II diabetic patients. The mathematical properties of the dynamic systems of glucose and insulin have been analyzed by some authors [2] [10] and [20] . The use of ordinary and partial differential equations to model biological systems cannot capture the rich variety of dynamics observed in natural systems for attempting to better our understanding of more and more complicated phenomena. One of the ways of dealing with these complexities is to include delays in the mathematical models. In 2012, a mathematical model has been developed to capture the integral impact of physical activity to glucose and insulin [21] . The existing global mathematical model due to John Thomas Sorensen [19] is complicated in term of computational complexity as it has many equations and parameters. The objective of this paper is therefore to propose a simplified mathematical model which is adaptable to observed data and can be used to adjust diabetes mellitus. Furthermore, we propose three compartment-mathematical model with only six odes by introducing a tissue compartment instead of a six compartment-mathematical model with 22 ordinary differential equations from [19] . The mathematical modeling follows a mathematical design where model equations are obtained by considering Adolf Fick's law [22] and Boyle-Mariotte's law [23] . For validation of the designed model, we use data obtained by solving the global mathematical model from [19] .

The rest of the paper is organised as follows. In Section 2, we set mathematical model equations. The Section 3 deals with qualitative study. Estimation of model parameters is presented in Section 4 while Section 5 focuses on concluding remarks.

2. Model Equations

The Sorensen model [19] deals with the behavior of different organs in a healthy human body by explicitly defining an individual compartment for each organ. The model also takes into account the individual effects of glucagon hormones and insulin on glucose metabolic rate in the different human body organ compartments. There are different methods that can be used to find mathematical model of glucose-insulin by physiologic compartmentalization. From the mathematical mode developed by Sorensen [19] , we take into account the work described in [24] and [25] . The current work modified the Sorensen model to consider the tissue compartment as set of brain, kidney, gut and periphery. This means that the blood from those elements passes through the tissue compartment to flow into the heart and lungs compartment. Figure 1 shows the exchange between physiologic compartments.

The metabolic sources and sinks in the glucose-insulin mathematical model are from the physiologic processes that happen at a constant rate. The mathematical nomenclature is defined in Table 1 and corresponds to the symbols used in the glucose insulin mathematical as shown in Figure 1. Since pancreatic insulin is released into the portal system which perfuses the liver, and since separate compartments have not been included in the model for vessel blood volumes, pancreatic insulin release appears as a source term in the liver insulin compartment.

Taking into account the exchanges illustrated in Figure 1, mass balances for the glucose insulin model result in a set of 6 simultaneous ordinary differential equations which are nonlinear as a result of the metabolic source and sink rates.

$\{\begin{array}{l}{V}_H^G\frac{\text{d}{G}_H\left(t\right)}{\text{d}t}=Q_L^G{G}_L\left(t\right)+{\gamma }_T^G{\left(G_T\left(t\right)\right)}^{\alpha }-Q_H^G{G}_H\left(t\right)-R_H^G,\hfill \\ V_L^G\frac{\text{d}{G}_L\left(t\right)}{\text{d}t}=Q_A^G{G}_H\left(t\right)-Q_L^G{G}_L\left(t\right)+R_L^G,\hfill \\ V_T^G\frac{\text{d}{G}_T\left(t\right)}{\text{d}t}=Q_P^G\left(G_H\left(t\right)-G_T\left(t\right)\right)-R_T^G,\hfill \\ V_H^I\frac{\text{d}{I}_H\left(t\right)}{\text{d}t}=Q_L^I{I}_L\left(t\right)+{\gamma }_T^I{\left(I_T\left(t\right)\right)}^{\beta }-Q_H^I{I}_H\left(t\right),\hfill \\ V_L^I\frac{\text{d}{I}_L\left(t\right)}{\text{d}t}=Q_A^I{I}_H\left(t\right)-Q_L^I{I}_L\left(t\right)+R_{PIR}-R_L^I,\hfill \\ V_T^I\frac{\text{d}{I}_T\left(t\right)}{\text{d}t}=Q_P^I\left(I_H\left(t\right)-I_T\left(t\right)\right)-R_T^I.\hfill \end{array}$ (1)

3. Model Qualitative Study

Let $X^e={\left(G_H^e,I_H^e,G_L^e,I_L^e,G_T^e,I_T^e\right)}^{\prime }$ be the steady state vector where X' denotes the transpose of X. In order to analyse the steady state, we need to solve the following system:

$\{\begin{array}{l}{Q}_L^G{G}_L^e+{\gamma }_T^G{\left(G_T^e\right)}^{\alpha }-Q_H^G{G}_H^e=R_H^G,\hfill \\ Q_A^G{G}_H^e-Q_L^G{G}_L^e=-R_L^G,\hfill \\ G_H^e-G_T^e=\frac{R_T^G}{Q_P^G},\hfill \\ Q_L^I{I}_L^e+{\gamma }_T^I{\left(I_T^e\right)}^{\beta }-Q_H^I{I}_H^e=0,\hfill \\ Q_A^I{I}_H^e-Q_L^I{I}_L^e=R_L^I-R_{PIR},\hfill \\ I_H^e-I_T^e=\frac{R_T^I}{Q_P^I}.\hfill \end{array}$ (2)

Note that the first three equations and the last three equations of (2) form the glucose model and insulin model respectively. From the glucose model we get $G_H^e$ and $G_L^e$ as functions of $G_T^e$

$G_H^e=G_T^e+\frac{R_T^G}{Q_P^G}\text{and}{G}_L^e=\frac{1}{Q_L^G}\left[-R_L^G+Q_A^G\left(G_T^e+\frac{R_T^G}{Q_P^G}\right)\right],$ (3)

and the glucose model becomes

$\{\begin{array}{l}{Q}_L^G{G}_L^e+{\gamma }_T^G{\left(G_T^e\right)}^{\alpha }-Q_H^G{G}_H^e=R_H^G,\hfill \\ G_H^e=G_T^e+\frac{R_T^G}{Q_P^G},\hfill \\ G_L^e=\frac{1}{Q_L^G}\left[-R_L^G+Q_A^G\left(G_T^e+\frac{R_T^G}{Q_P^G}\right)\right].\hfill \end{array}$

In the same way, from the insulin model we get $I_H^e$ and $I_L^e$ as functions of $I_T^e$

$I_H^e=I_T^e+\frac{R_T^I}{Q_P^I}\text{and}{I}_L^e=\frac{1}{Q_L^I}\left[Q_A^I\left(I_T^e+\frac{R_T^I}{Q_P^I}\right)+R_{PIR}-R_L^I\right],$ (4)

and the insulin model can be rewritten as follows

$\{\begin{array}{l}{Q}_L^I{I}_L^e+{\gamma }_T^I{\left(I_T^e\right)}^{\beta }-Q_H^I{I}_H^e=0,\hfill \\ I_H^e=I_T^e+\frac{R_T^I}{Q_P^I},\hfill \\ I_L^e=\frac{1}{Q_L^I}\left[Q_A^I\left(I_T^e+\frac{R_T^I}{Q_P^I}\right)+R_{PIR}-R_L^I\right].\hfill \end{array}$

Let ${\mathcal{J}}_G=\frac{\partial G}{\partial X}\left(X^e\right)$ and ${\mathcal{J}}_I=\frac{\partial I}{\partial X}\left(X^e\right)$ be Jacobian matrices of glucose model and insulin model respectively where all derivatives are evaluated at the equilibrium point $X^e$ . After some algebraic calculations we get

${\mathcal{J}}_G=\left(\begin{array}{ccc}-Q_H^G& Q_L^G& \alpha {\gamma }_T^G{\left(G_T^e\right)}^{\alpha -1}\\ Q_A^G& -Q_L^G& 0\\ Q_P^G& 0& -Q_P^G\end{array}\right)\text{and}{\mathcal{J}}_I=\left(\begin{array}{ccc}-Q_H^I& Q_L^I& \beta {\gamma }_T^I{\left(I_T^e\right)}^{\beta -1}\\ Q_A^I& -Q_L^I& 0\\ Q_P^I& 0& -Q_P^I\end{array}\right)\mathrm{.}$ (5)

The behaviour of the mathematical model (1) near the steady state can be analysed by the nature of the real parts of the eigenvalues of matrices ${\mathcal{J}}_G$ and ${\mathcal{J}}_I$ . The proof of the theorem below will use the following proposition due to Routh-Hurwitz [26] .

Proposition 1. Let $a_1\mathrm{,}{a}_2$ and $a_3$ be positive real numbers. The roots of the polynomial

${\lambda }^3+a_1{\lambda }^2+a_2\lambda +a_3=0$

have negative real parts when $a_1{a}_2>a_3$ .

Theorem 2.

The system (1) is asymptotically stable if the following conditions are satisfied:

1)

$\alpha {\gamma }_T^G{\left(G_T^e\right)}^{\alpha -1}<Q_H^G-Q_A^G,$ (6)

2)

$\alpha {\gamma }_T^G{\left(G_T^e\right)}^{\alpha -1}<\frac{Q_L^G}{Q_P^G}\left(Q_H^G-Q_A^G\right)+Q_H^G+Q_L^G,$ (7)

3)

$\begin{array}{c}\alpha {\gamma }_T^G{\left(G_T^e\right)}^{\alpha -1}<\frac{1}{\left[Q_P^G\left(Q_P^G+Q_H^G\right)\right]}[2Q_P^G{Q}_L^G{Q}_H^G-Q_L^G{Q}_A^G{Q}_H^G+{\left(Q_H^G\right)}^2\left(Q_P^G+Q_L^G\right)\\ +{\left(Q_L^G\right)}^2\left(Q_P^G+Q_H^G-Q_A^G\right)+{\left(Q_P^G\right)}^2\left(Q_H^G+Q_L^G\right)],\end{array}$ (8)

4)

$\beta {\gamma }_T^I{\left(I_T^e\right)}^{\beta -1}<Q_H^I-Q_A^I,$ (9)

5)

$\beta {\gamma }_T^I{\left(I_T^e\right)}^{\beta -1}<\frac{Q_L^I}{Q_P^I}\left(Q_H^I-Q_A^I\right)+Q_H^I+Q_L^I,$ (10)

6)

$\begin{array}{c}\beta {\gamma }_T^I{\left(I_T^e\right)}^{\beta -1}<\frac{1}{\left[Q_P^I\left(Q_P^I+Q_H^I\right)\right]}[2Q_P^I{Q}_L^I{Q}_H^I-Q_L^I{Q}_A^I{Q}_H^I+{\left(Q_H^I\right)}^2\left(Q_P^I+Q_L^I\right)\\ +{\left(Q_L^I\right)}^2\left(Q_P^I+Q_H^I-Q_A^I\right)+{\left(Q_P^I\right)}^2\left(Q_H^I+Q_L^I\right)].\end{array}$ (11)

Proof.

The system (1) is asymptotically stable if and only if ${\mathcal{J}}_G$ and ${\mathcal{J}}_I$ are stability matrix; that is, every eigenvalue of ${\mathcal{J}}_G$ and ${\mathcal{J}}_I$ has a negative real part. The characteristic equation associated to ${\mathcal{J}}_G$ is $P_G\left(\lambda \right)=0$ given by

$\left|\begin{array}{ccc}-Q_H^G-\lambda & Q_L^G& \alpha {\gamma }_T^G{\left(G_T^e\right)}^{\alpha -1}\\ Q_A^G& -Q_L^G-\lambda & 0\\ Q_P^G& 0& -Q_P^G-\lambda \end{array}\right|=0.$

That is

$\begin{array}{l}{\lambda }^3+\left(Q_H^G+Q_L^G+Q_P^G\right){\lambda }^2\\ +\left[Q_P^G\left(Q_H^G+Q_L^G-\alpha {\gamma }_T^G{\left(G_T^e\right)}^{\alpha -1}\right)+Q_L^G\left(Q_H^G-Q_A^G\right)\right]\lambda \\ +Q_P^G{Q}_L^G\left(Q_H^G-Q_A^G-\alpha {\gamma }_T^G{\left(G_T^e\right)}^{\alpha -1}\right)=0.\end{array}$ (12)

Similarly $P_I\left(\lambda \right)=0$ is equivalent to

$\begin{array}{l}{\lambda }^3+\left(Q_H^I+Q_L^I+Q_P^I\right){\lambda }^2\\ +\left[Q_P^I\left(Q_H^I+Q_L^I-\alpha {\gamma }_T^I{\left(I_T^e\right)}^{\alpha -1}\right)+Q_L^I\left(Q_H^I-Q_A^I\right)\right]\lambda \\ +Q_P^I{Q}_L^I\left(Q_H^I-Q_A^I-\alpha {\gamma }_T^I{\left(I_T^e\right)}^{\alpha -1}\right)=0.\end{array}$ (13)

Using Proposition 1, we need to verify the following requirements: all the coefficients of (12) are positive and the product of coefficients of second and third terms of (12) is strictly greater than its fourth term.

Indeed, since all vascular blood flow rates are positive, then the relation $Q_H^G+Q_L^G+Q_P^G>0$ is obvious. The next requirement is

$Q_P^G{Q}_L^G\left(Q_H^G-Q_A^G-\alpha {\gamma }_T^G{\left(G_T^e\right)}^{\alpha -1}\right)>0,$

which after some calculations is equivalent to

$\alpha {\gamma }_T^G{\left(G_T^e\right)}^{\alpha -1}<Q_H^G-Q_A^G.$

Similarly, the requirement

$Q_P^G\left(Q_H^G+Q_L^G-\alpha {\gamma }_T^G{\left(G_T^e\right)}^{\alpha -1}\right)+Q_L^G\left(Q_H^G-Q_A^G\right)>0,$

is equivalent to

$\alpha {\gamma }_T^G{\left(G_T^e\right)}^{\alpha -1}<\frac{Q_L^G}{Q_P^G}\left(Q_H^G-Q_A^G\right)+Q_H^G+Q_L^G.$

The last requirement

$\begin{array}{l}\left[Q_P^G\left(Q_H^G+Q_L^G-\alpha {\gamma }_T^G{\left(G_T^e\right)}^{\alpha -1}\right)+Q_L^G\left(Q_H^G-Q_A^G\right)\right]\left(Q_H^G+Q_L^G+Q_P^G\right)\\ >Q_P^G{Q}_L^G\left(Q_H^G-Q_A^G-\alpha {\gamma }_T^G{\left(G_T^e\right)}^{\alpha -1}\right)\mathrm{,}\end{array}$ (14)

yields after calculations (8).

The requirements (9), (10) and (11) are obtained in a similar way by considering the insulin model.,

4. Model Parameter Identification

The nonlinear system (1) can be represented in the following compact form

$\stackrel{\dot{}}{X}\left(t\right)=f\left(X\left(t\right),\mu \right),$ (15)

where $\mu$ is the vector of parameters to be estimated. That is

$\mu ={\left({\gamma }_T^G,\alpha ,R_H^G,\beta ,{\gamma }_T^I,R_L^G,R_T^G,R_T^I\right)}^{\prime }.$

The mathematical model requires parameter identification which can be carried out by setting the following optimal control problem. We determine the control $\mu$ such that the cost functional

$J\left(\mu \right)={\displaystyle {\int }_0^{t_f}}{q}_G{\left(G_H\left(t\right)-G_H^{obs}\right)}^2+q_I{\left(I_H\left(t\right)-I_H^{obs}\right)}^2\text{d}t,$ (16)

is minimized under the restriction of the model Equation (15). The positive scalar coefficients $q_G$ and $q_I$ determine how much weight is associated to each term in the integrand. Superscript “obs” refers to the observed state to which the system is transferred by the control. In order to obtain the observed data, we solve the global model of [19] . The variation of glucose and insulin in heart are plotted in Figure 2.

For computational purposes we discretize the system (15) using N linear B-splines. Let us consider

${\mathcal{B}}^N=\left\{{\psi }_j^N,j=1,\cdots ,N\right\},$ (17)

a linear B-splines basis functions on the uniform grid

${\Omega }_N=\left\{t_k=\frac{kT}{N},k=0,\cdots ,N\right\},$ (18)

such that

${\psi }_i^N\left(t_k\right)={\delta }_{ik},$

where ${\delta }_{ik}$ is Kronecker symbol. Let us introduce the vector space $W^N$ whose

the basis is ${\mathcal{B}}^N$ . It follows that $\dim W^N=N$ and $W^n\subset W^{n+1},n=1,\cdots ,N$ . We assume that functions appearing in the system (15) are continuous on $\left[\mathrm{0,}T\right]$ . Let us denote $W=C^0\left(0,T\right)$ and consider the interpolation operator

${\Pi }^N:W\to W^N\mathrm{,}$

satisfying $\forall \varphi \in W$

${\Pi }^N\varphi \left(t_k\right)=\varphi \left(t_k\right),\mathrm{<mi>\; </mi>\; <mi>\; </mi>\; <mi>\; </mi>}k=1,\cdots ,N.$

Therefore, in this setting we are looking for the solution $X^N\in W^N$ of the following discrete problem

${\stackrel{\dot{}}{X}}^N\left(t\right)=f\left(X^N\left(t\right)\mathrm{,}{\mu }^N\right)\mathrm{,}\text{such}\text{that}{X}^N\left(0\right)=X_0^N\mathrm{,}$ (19)

where the control is

${\mu }^N=\left({\gamma }_T^{G,N},{\alpha }^N,R_H^{G,N},{\beta }^N,{\gamma }_T^{I,N},R_L^{G,N},R_T^{G,N},R_T^{I,N}\right)\in {\left(W^N\right)}^8.$

The corresponding discrete optimal control problem (16) is to minimize

${\displaystyle \underset{k=1}{\overset{N}{\sum }}}\left(q_G{\left(G_H^N\left(t_k\right)-G_H^{obs}\right)}^2+q_I{\left(I_H^N\left(t_k\right)-I_H^{obs}\right)}^2\right)h\text{with}h=\frac{T}{N}$ (20)

with respect to (19). In compact form the problem (20) can be rewritten as follows

$\underset{{\mu }^N}{\min }{J}^N\left({\mu }^N\right)={\displaystyle \underset{k=1}{\overset{N}{\sum }}}h{Y}_k^{\text{T}}R{Y}_k$ (21)

subject to

$\{\begin{array}{l}{\stackrel{\dot{}}{X}}^N\left(t\right)=f\left(X^N\left(t\right),{\mu }^N\right)\hfill \\ X^N\left(0\right)=X_0^N\hfill \end{array}$ (22)

where $Y_k,k=1,\cdots ,N$ is the following matrix

${\left(\left(G_H^N\left(t_k\right)-G_H^{obs}\right),\left(I_H^N\left(t_k\right)-I_H^{obs}\right)\right)}^{\prime },$

and R is the matrix defined by

$R=\left(\begin{array}{cc}{q}_G& 0\\ 0& q_I\end{array}\right).$

The numerical computations have been carried out using a collection of MaTlaB routines [27] specifically the built-in function fmincon. Table 2 shows the estimated parameters of mathematical model (1). Using those parameters we get numerical solutions shown in Figure 3.

5. Concluding Remarks

Physiological and dynamical conditions for glucose and insulin impose the need for relatively simple models that should be able to describe as accurately as possible the mechanical behavior of glucose-insulin system. In this work we have proposed a three compartmental mathematical model that describes the variation of glucose and insulin for human being. The modelling technique is used to

provide interesting answers to the question of determining the global mathematical model with lower number of equations for glucose-insulin system. Numerical results show that the proposed model is adaptable to data. In fact Figure 3 shows that adjustment for glucose and insulin trends match observed data. The proposed mathematical model can also be used to adjust diabetes mellitus type I or type II for diabetic patients.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

Cite this paper

Ntaganda, J.M., Minani, F., Banzi, W., Mpinganzima, L., Niyobuhungiro, J., Gahutu, J.B., Rutaganda, E., Kambutse, I. and Dusabejambo, V. (2018) Simplified Mathematical Model of Glucose-Insulin System. American Journal of Computational Mathematics, 8, 233-244. https://doi.org/10.4236/ajcm.2018.83019

References