Dynamics of a glucose–insulin model

Abstract

Diabetes mellitus is a noncommunicable disease, which is a serious threat to human health around the world. In this paper, we propose a simple glucose–insulin model with Michaelis–Menten function as insulin degradation rate to mimic the pathogenic mechanism of diabetes. By theoretical analysis, a unique positive equilibrium of model exists and it is globally asymptotically stable. The four strategies are designed for diabetes patients based on the sensitivity of parameters, including insulin injection and medicine treatments. Numerical simulations are given to support the theoretical results.

Keywords:

• Diabetes
• glucose
• insulin
• therapy
• equilibrium

1. Introduction

Diabetes mellitus is a chronic and complex metabolic illness characterized by high blood glucose level for long duration, which is one of the prevalent diseases affecting millions of persons around the world nowadays. Long-term hyperglycemia also leads to many other problems of organs of the body, such as cardiovascular disease, diabetes feet, cataract, nephropathy and neuropathy and so on. By the World Health Organization (WHO) and the International Diabetes Federation (IDF) criterion, there are four types: Type 1 diabetes (T1DM), Type 2 diabetes (T2DM), gestational diabetes and other specific types of diabetes, in which the T1DM and T2DM account for above 95%. T1DM is considered to be an autoimmune disorder because the β cells in the pancreas are damaged and cannot secrete insulin. T2DM is considered to be the result of insulin resistance, which the β cells cannot produce enough insulin relatively or the secreted insulin has no ability to decompose glucose correctly [1]. In view of large proportion of T1DM and T2DM, we only focus on the two types in the paper.

How to judge a person to get the diabetes mellitus? Generally speaking, early typical symptoms include polyuria, thirst, more food and weight loss. The common glucose diagnostic tests in clinical practice are intravenous glucose tolerance test (IVGTT) and oral glucose tolerance test (OGTT). IVGTT is a test in which glucose is injected intravenously and blood samples are collected following the glucose injection [2], and OGTT is a test in which glucose is given orally and blood samples are collected over 2 hours following the glucose infusion [3]. If fasting blood glucose value is greater than 7 $\mathrm{mmol}/\mathrm{l}$ and the value is higher than 11.1 $\mathrm{mmol}/\mathrm{l}$ 2 hours after meals, it can be identified as diabetes mellitus. From these tests, doctors can obtain some indexes, such as insulin sensitivity, glucose effectiveness, insulin action, β cells function and so on, which can be used to diagnose disease and give some treatment measures. The current therapies are mainly hypodermic insulin injections and oral hypoglycemic drugs, which could be used to treat the symptoms, but the character of diabetes dynamics is still not fully understood. Therefore, the work of studying diabetes is necessary and urgent.

In the last several decades, many mathematical models about glucose-insulin regulatory system have been developed and analysed, and some good results have been obtained [4–10]. Bolie [4] is considered as the pioneer to study glucose–insulin dynamics. He developed a minimal model to mimic blood glucose regulation. The model was written as $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}x}{\mathrm{d}t}=p-\alpha x+\beta y}\\ {\displaystyle \frac{\mathrm{d}y}{\mathrm{d}t}=q-\gamma x-\delta y,}\end{array}$where x and y represent the deviation in insulin and glucose concentration from mean physiological value, α denotes the sensitivity of insulinase activity to elevate insulin concentration, β denotes the sensitivity of pancreatic insulin to elevate glucose concentration, δ and γ represent the combined sensitivity of liver glycogen storage and tissue glucose utilization to elevate insulin and glucose concentration respectively. The author studied the model and evaluated the values of four coefficients and their biological variations, and the results were used to manage the diabetes.

Based on Bolie's work, a lot of different kinds of mathematical models have appeared for studying the glucose–insulin regulatory system. Gaetano and Arino [5] proposed a ‘dynamical model’ with delay. They studied the stability of unique equilibrium and gave the positive and bounded of solutions. The dynamical model overcame the limitations of coupled minimal model and was more attractive physiologically because of delay. [6–10] also studied the glucose–insulin models with delay and obtained some good conclusions. Lombarte et al. [11] developed a glucose–insulin model for the homeostasis of glucose in diabetic rats. They obtained corresponding parameters from the values of glucose and insulin in blood and urine. For more related studies, see [12–14] and references therein.

Lots of good progress have been made in studying the glucose–insulin regulatory systems, but there are still some shortcomings. For example, some parameters are not identifiable or fit to experimental data poorly [15]. Generally speaking, once diagnosed with diabetes, the patients will be advised to inject the insulin for controlling the level of glucose. However, how doctors decide the level of insulin injection? There is no uniform available formula. To solve this problem, many works were done for searching suitable strategies. For example, Maianti et al. [16] hypothesized using insulin degradation enzyme to improve glucose tolerance for T2DM. Huard et al. [17] studied a deterministic nonlinear model with two delays and explored four strategies for restoring healthy regulation. In this paper, we expect to find some methods to maintain glucose in normal level by mathematical models and obtain more efficient and economic therapies for diabetics.

This paper is organized as follows: in Section 2, the model about glucose and insulin will be proposed. In Section 3, the basic properties of the model will be given, such as positivity and boundedness of solutions. The conditions of local and global stability of equilibrium are established applying characteristic roots and Lyapunov function in Section 4. The sensitivity of parameters is studied in Section 5 and four corresponding strategies are proposed in Section 6. In the last section, we discuss our results and make a brief summary.

2. Model

As well known, the typical exogenous infusions of glucose come from carbohydrate in daily life, which include cereals and meat ingestion and so on. They afford energy to maintain all kinds of complex life activities. The consumption ways of glucose have two aspects: one is insulin-independent glucose consumers, such as brain, another is insulin-dependent glucose consumers, such as muscle and so on. The insulin includes synthetic insulin and secreted by β cells in pancreas. The raised blood glucose level stimulates β cells in the pancreas to secret insulin which helps the decomposition of glucose. In this paper, we use $G(t)$ denotes plasma glucose concentration (mg/dl) and $I(t)$ denotes plasma insulin concentration (μU/ml). $G_{in}$ and $I_{in}$ denote the glucose intake rate and the exogenous insulin infusion rate which are all constants, respectively.

In existing glucose–insulin models, the response functions of insulin clearance are assumed to be proportional to insulin concentration [7–10]. The assumption may be inadequate for looking for more appropriate diagnosis and cure strategies. As pointed out in [10], the change rate increases with the increasing resource but does not increase indefinitely even if the resource is excess. So they considered an enhanced modelling of the glucose–insulin system, in which the insulin degradation rate is Michaelis–Menten function $\frac{dI}{e+I}$, where d is the maximum insulin clearance rate and e is the half-saturation value. It is more realistic than linear rate. Motivated by the work of [10], a glucose–insulin model is built. The model structure is shown in Figure 1, and the complete model is given as follows. (1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}G}{\mathrm{d}t}=G_{in}-aG-bGI,}\\ {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=I_{in}+cG-\frac{dI}{e+I},}\end{array}$(1) with initial conditions $G(0)>0,I(0)>0$. The meaning of parameters is showed in Table 1.

Table 1. The state variables for model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}G}{\mathrm{d}t}=G_{in}-aG-bGI,}\\ {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=I_{in}+cG-\frac{dI}{e+I},}\end{array}$(1) ).

Parameters	Description	Unit
$G_{in}$	The glucose intake rate	mg/dl/min
$I_{in}$	The exogenous insulin infusion rate	μU/ml/min
a	The insulin-independent utilization rate	min${}^{-1}$
b	The insulin-dependent utilization rate(insulin resistance)	ml/μU/min
c	The insulin secretion rate	μU/ml/min/(mg/dl)
d	The maximum insulin clearance rate	μU/ml/min
e	The half-saturation value	–

Figure 1. The glucose–insulin system diagram.

3. The basic properties of model (1)

In view of the biological meaning, the variables of the model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}G}{\mathrm{d}t}=G_{in}-aG-bGI,}\\ {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=I_{in}+cG-\frac{dI}{e+I},}\end{array}$(1) ) must be positive for all time. In the following, the positivity and boundedness of solutions of model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}G}{\mathrm{d}t}=G_{in}-aG-bGI,}\\ {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=I_{in}+cG-\frac{dI}{e+I},}\end{array}$(1) ) will be given.

To study the properties of the steady state of model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}G}{\mathrm{d}t}=G_{in}-aG-bGI,}\\ {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=I_{in}+cG-\frac{dI}{e+I},}\end{array}$(1) ), the equilibria of system should be determined first. Let the right of model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}G}{\mathrm{d}t}=G_{in}-aG-bGI,}\\ {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=I_{in}+cG-\frac{dI}{e+I},}\end{array}$(1) ) is equal to zero, that is (2) $\{\begin{array}{l}{G}_{in}-aG-bGI=0,\\ {\displaystyle I_{in}+cG-\frac{dI}{e+I}=0.}\end{array}$(2) When $d>I_{in}$, there exists a unique positive equilibrium $(G^{\ast },I^{\ast })$, where $G^{\ast }=\frac{G_{in}}{a+b{I}^{\ast }},I^{\ast }=\frac{(ad-c{G}_{in}-a{I}_{in}-be{I}_{in})-\sqrt{\mathrm{\Delta }}}{2b(I_{in}-d)}$and $\mathrm{\Delta }=(ad-c{G}_{in}-a{I}_{in}-be{I}_{in}{)}^2-4be(I_{in}-d)(c{G}_{in}+a{I}_{in}).$

Lemma 3.1

For model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}G}{\mathrm{d}t}=G_{in}-aG-bGI,}\\ {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=I_{in}+cG-\frac{dI}{e+I},}\end{array}$(1) ), all solutions with initial conditions are positive; furthermore, all solutions are bounded for t>0 if $G(0)\le \frac{G_{in}}{a}$ and $d>I_{in}+\frac{c{G}_{in}}{a}$.

Proof.

Let $(G,I)$ is the solution of model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}G}{\mathrm{d}t}=G_{in}-aG-bGI,}\\ {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=I_{in}+cG-\frac{dI}{e+I},}\end{array}$(1) ). If $G(t)\le 0$, then there exists a $t_0>0$ such that $G(t_0)=0$ and $G^{\prime }(t_0)<0$.

But ${\frac{\mathrm{d}G}{\mathrm{d}t}|}_{t=t_0}=G_{in}-aG(t_0)-bG(t_0)I(t_0)=G_{in}>0,$this is a contradiction. So $G(t)>0$ for t>0. By the same argument, $I(t)>0$ for t>0.

From model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}G}{\mathrm{d}t}=G_{in}-aG-bGI,}\\ {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=I_{in}+cG-\frac{dI}{e+I},}\end{array}$(1) ), we have $\frac{\mathrm{d}G}{\mathrm{d}t}=G_{in}-aG-bGI\le G_{in}-aG,$so $G\le \frac{G_{in}}{a}$ if $G(0)\le \frac{G_{in}}{a}$, that is $G(t)$ is bounded for all t>0.

The second equation of model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}G}{\mathrm{d}t}=G_{in}-aG-bGI,}\\ {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=I_{in}+cG-\frac{dI}{e+I},}\end{array}$(1) ) is changed as $\frac{\mathrm{d}I}{\mathrm{d}t}=I_{in}+cG-\frac{dI}{e+I}\le I_{in}+\frac{c{G}_{in}}{a}-\frac{dI}{e+I}.$Obviously, if $I>\frac{e(I_{in}+\frac{c{G}_{in}}{a})}{d-(I_{in}+\frac{c{G}_{in}}{a})}$, then $I^{\prime }(t)<0$. In fact, by Lemma 2.2 in [10], with the condition $d>I_{in}+\frac{c{G}_{in}}{a}$, we can get $\underset{t\to \mathrm{\infty }}{lim sup}I(t)\le \frac{e(I_{in}+\frac{c{G}_{in}}{a})}{d-(I_{in}+\frac{c{G}_{in}}{a})},$therefore $I(t)$ is bounded for all t>0. The proof is complete.

4. The stability of equilibrium

In this section, the local and global stabilities of equilibrium of model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}G}{\mathrm{d}t}=G_{in}-aG-bGI,}\\ {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=I_{in}+cG-\frac{dI}{e+I},}\end{array}$(1) ) will be established.

Suppose $(G,I)$ is an arbitrary equilibrium of model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}G}{\mathrm{d}t}=G_{in}-aG-bGI,}\\ {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=I_{in}+cG-\frac{dI}{e+I},}\end{array}$(1) ), then the Jacobian matrix of model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}G}{\mathrm{d}t}=G_{in}-aG-bGI,}\\ {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=I_{in}+cG-\frac{dI}{e+I},}\end{array}$(1) ) at $(G,I)$ is $D=\left(\begin{array}{cc}-a-bI& -bG\\ c& {\displaystyle -\frac{de}{(e+I{)}^2}}\end{array}\right),$the corresponding characteristic equation of the model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}G}{\mathrm{d}t}=G_{in}-aG-bGI,}\\ {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=I_{in}+cG-\frac{dI}{e+I},}\end{array}$(1) ) at $(G^{\ast },I^{\ast })$ is $\begin{array}{rl}|\lambda E-D{|}_{(G^{\ast },I^{\ast })}& =(\lambda +a+b{I}^{\ast })(\lambda +\frac{de}{(e+I^{\ast }{)}^2})+cb{G}^{\ast }\\ & ={\lambda }^2+(a+b{I}^{\ast }+\frac{de}{(e+I^{\ast }{)}^2})\lambda +(a+b{I}^{\ast })\frac{de}{(e+I^{\ast }{)}^2}+cb{G}^{\ast }\\ & =0,\end{array}$then ${\lambda }_{1,2}=\frac{-(a+b{I}^{\ast }+\frac{de}{(e+I^{\ast }{)}^2})\pm \sqrt{(a+b{I}^{\ast }+\frac{de}{(e+I^{\ast }{)}^2}{)}^2-4((a+b{I}^{\ast })\frac{de}{(e+I^{\ast }{)}^2}+bc{G}^{\ast })}}{2}.$It is easy to obtain that two characteristic roots have negative real parts, then the unique positive equilibrium $(G^{\ast },I^{\ast })$ is locally asymptotically stable if $d>I_{in}$.

Next, we will establish the global asymptotical stability of equilibrium $(G^{\ast },I^{\ast })$ under the existence condition by constructing a suitable Lyapunov function.

Theorem 4.1

The positive equilibrium $(G^{\ast },I^{\ast })$ of model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}G}{\mathrm{d}t}=G_{in}-aG-bGI,}\\ {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=I_{in}+cG-\frac{dI}{e+I},}\end{array}$(1) ) is globally asymptotically stable if $d>I_{in}$.

Proof.

Define a Lyapunov function as follows: $L(t)=\frac{c}{2}(G-G^{\ast }{)}^2+\frac{b{G}^{\ast }}{2}(I-I^{\ast }{)}^2,$then we calculate the derivative of $L(t)$ along the solution of model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}G}{\mathrm{d}t}=G_{in}-aG-bGI,}\\ {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=I_{in}+cG-\frac{dI}{e+I},}\end{array}$(1) ). $\begin{array}{rl}\dot{L}(t)& =c(G-G^{\ast })\dot{G}+b{G}^{\ast }(I-I^{\ast })\dot{I}\\ & =c(G-G^{\ast })(G_{in}-aG-bGI)+b{G}^{\ast }(I-I^{\ast })(I_{in}+cG-\frac{dI}{e+I})\\ & =c(G-G^{\ast })(a(G^{\ast }-G)+b(G^{\ast }{I}^{\ast }-GI)\\ & +b{G}^{\ast }(I-I^{\ast })(c(G-G^{\ast })-\frac{de(I^{\ast }-I)}{(e+I)(e+I^{\ast })})\\ & =-ac(G-G^{\ast }{)}^2-bcI(G-G^{\ast }{)}^2-\frac{bde{G}^{\ast }(I-I^{\ast }{)}^2}{(e+I)(e+I^{\ast })}\\ & \le 0.\end{array}$It is obvious that $\dot{L}(t)=0$ holds if and only if $G=G^{\ast }$ and $I=I^{\ast }$, then the unique positive equilibrium $(G^{\ast },I^{\ast })$ is globally asymptotically stable if $d>I_{in}$.

5. The sensitivity of parameters

Sensitivity of parameters is a useful tool to study the variations of parameters qualitatively and quantitatively. To find out the effect of parameters on the equilibrium $(G^{\ast },I^{\ast })$, we first discuss the relationships of parameters $(b,c,d)$ and equilibrium $(G^{\ast },I^{\ast })$. The method is differentiating implicitly Equation (2(2) $\{\begin{array}{l}{G}_{in}-aG-bGI=0,\\ {\displaystyle I_{in}+cG-\frac{dI}{e+I}=0.}\end{array}$(2) ) with respect to $b,c$ and d, respectively. We always denote the derivative x with respect to t as $x_t$.

5.1. The relationships of b, c, d and $(G^{\ast },I^{\ast })$

$\{\begin{array}{l}-a{G}_b^{\ast }-G^{\ast }{I}^{\ast }-b(G_b^{\ast }{I}^{\ast }+G^{\ast }{I}_b^{\ast })=0,\\ {\displaystyle c{G}_b^{\ast }-\frac{de{I}_b^{\ast }}{(e+I^{\ast }{)}^2}=0,}\end{array}$that is $\left(\begin{array}{cc}c& {\displaystyle -\frac{de}{(e+I^{\ast }{)}^2}}\\ a+b{I}^{\ast }& b{G}^{\ast }\end{array}\right)\left(\begin{array}{c}{G}_b^{\ast }\\ I_b^{\ast }\end{array}\right)=\left(\begin{array}{c}0\\ -G^{\ast }{I}^{\ast }\end{array}\right).$Since ${\mathrm{\Delta }}_1=bc{G}^{\ast }+\frac{de}{(e+I^{\ast }{)}^2}(a+b{I}^{\ast })>0$, then $\left(\begin{array}{c}{G}_b^{\ast }\\ I_b^{\ast }\end{array}\right)=\frac{1}{{\mathrm{\Delta }}_1}\left(\begin{array}{cc}b{G}^{\ast }& {\displaystyle \frac{de}{(e+I^{\ast }{)}^2}}\\ -(a+b{I}^{\ast })& c\end{array}\right)\left(\begin{array}{c}0\\ -G^{\ast }{I}^{\ast }\end{array}\right)=\frac{1}{{\mathrm{\Delta }}_1}\left(\begin{array}{c}{\displaystyle -\frac{de{G}^{\ast }{I}^{\ast }}{(e+I^{\ast }{)}^2}}\\ -c{G}^{\ast }{I}^{\ast }\end{array}\right).$The result indicates $G^{\ast }$ and $I^{\ast }$ will decrease with increased b, which implies that improving the insulin-dependent utilization rate could reduce the glucose $G^{\ast }$ and insulin $I^{\ast }$. $\{\begin{array}{l}a{G}_c^{\ast }+b{G}_c^{\ast }{I}^{\ast }+b{G}^{\ast }{I}_c^{\ast }=0,\\ {\displaystyle G^{\ast }+c{G}_c^{\ast }-\frac{de{I}_c^{\ast }}{(e+I^{\ast }{)}^2}=0,}\end{array}$that is $\left(\begin{array}{cc}a+b{I}^{\ast }& b{G}^{\ast }\\ c& {\displaystyle -\frac{de}{(e+I^{\ast }{)}^2}}\end{array}\right)\left(\begin{array}{c}{G}_c^{\ast }\\ I_c^{\ast }\end{array}\right)=\left(\begin{array}{c}0\\ -G^{\ast }\end{array}\right).$Since ${\mathrm{\Delta }}_2=-(a+b{I}^{\ast })\frac{de}{(e+I^{\ast }{)}^2}-bc{G}^{\ast }=-{\mathrm{\Delta }}_1<0$, then $\left(\begin{array}{c}{G}_c^{\ast }\\ I_c^{\ast }\end{array}\right)=\frac{1}{{\mathrm{\Delta }}_1}\left(\begin{array}{cc}{\displaystyle -\frac{de}{(e+I^{\ast }{)}^2}}& -b{G}^{\ast }\\ -c& a+b{I}^{\ast }\end{array}\right)\left(\begin{array}{c}0\\ -G^{\ast }\end{array}\right)=-\frac{1}{{\mathrm{\Delta }}_1}\left(\begin{array}{c}b{G^{\ast }}^2\\ -G^{\ast }(a+b{I}^{\ast })\end{array}\right).$It is obvious that $G^{\ast }$ decreases with increased c, and $I^{\ast }$ increases with increased c. The result indicates if the secretion ability of insulin enhances, then the insulin level must increase and the glucose level becomes less. $\{\begin{array}{l}-a{G}_d^{\ast }-b(G_d^{\ast }{I}^{\ast }+G^{\ast }{I}_d^{\ast })=0,\\ {\displaystyle c{G}_d^{\ast }-\frac{e{I}^{\ast }+(I^{\ast }{)}^2+de{I}_d^{\ast }}{(e+I^{\ast }{)}^2}=0,}\end{array}$that is $\left(\begin{array}{cc}a+b{I}^{\ast }& b{G}^{\ast }\\ c(e+I^{\ast }{)}^2& -de\end{array}\right)\left(\begin{array}{c}{G}_d^{\ast }\\ I_d^{\ast }\end{array}\right)=\left(\begin{array}{c}0\\ I^{\ast }(e+I^{\ast })\end{array}\right).$Since ${\mathrm{\Delta }}_3=-de(a+b{I}^{\ast })-bc{G}^{\ast }(e+I^{\ast }{)}^2=-(e+I^{\ast }{)}^2{\mathrm{\Delta }}_1<0$, then $\begin{array}{rl}\left(\begin{array}{c}{G}_d^{\ast }\\ I_d^{\ast }\end{array}\right)& =-\frac{1}{(e+I^{\ast }{)}^2{\mathrm{\Delta }}_1}\left(\begin{array}{cc}-de& -b{G}^{\ast }\\ -c(e+I^{\ast }{)}^2& a+b{I}^{\ast }\end{array}\right)\left(\begin{array}{c}0\\ I^{\ast }(e+I^{\ast })\end{array}\right)\\ & =-\frac{1}{(e+I^{\ast }){\mathrm{\Delta }}_1}\left(\begin{array}{c}-b{G}^{\ast }{I}^{\ast }\\ (a+b{I}^{\ast })I^{\ast }\end{array}\right).\end{array}$According to the direct calculation, the level of $G^{\ast }$ will increase with increased d, but $I^{\ast }$ will increase with decreased d. Therefore, increasing the maximum insulin clearance rate can improve the glucose concentration $G^{\ast }$ and reduce the insulin concentration $I^{\ast }$.

5.2. Sensitivity indices of $(G^{\ast },I^{\ast })$

In the subsection, we discuss the sensitivity indices of $(G^{\ast },I^{\ast })$ to parameters $(b,c,d)$ by the method of [18,19]. The sensitivity indices can be used to measure the relative change of a state variable as a parameter changes.

Definition 5.1

The normalized forward sensitivity index of a variable u that depends differentiable on a parameter p is defined as ${\mathrm{{\rm Y}}}_p^u:=\frac{\mathrm{\partial }u}{\mathrm{\partial }p}\times \frac{p}{u}.$

Now, we calculate the normalized forward sensitivity indices as follows: $\begin{array}{rl}{\mathrm{{\rm Y}}}_b^{G^{\ast }}& =\frac{b}{G^{\ast }}\frac{\mathrm{\partial }{G}^{\ast }}{\mathrm{\partial }b}=-\frac{bde{I}^{\ast }}{{\mathrm{\Delta }}_1(e+I^{\ast }{)}^2},{\mathrm{{\rm Y}}}_b^{I^{\ast }}=\frac{b}{I^{\ast }}\frac{\mathrm{\partial }{I}^{\ast }}{\mathrm{\partial }b}=-\frac{bc{G}^{\ast }}{{\mathrm{\Delta }}_1},\\ {\mathrm{{\rm Y}}}_c^{G^{\ast }}& =\frac{c}{G^{\ast }}\frac{\mathrm{\partial }{G}^{\ast }}{\mathrm{\partial }c}=-\frac{bc{G}^{\ast }}{{\mathrm{\Delta }}_1},{\mathrm{{\rm Y}}}_c^{I^{\ast }}=\frac{c}{I^{\ast }}\frac{\mathrm{\partial }{I}^{\ast }}{\mathrm{\partial }c}=\frac{c(a+b{I}^{\ast })G^{\ast }}{{\mathrm{\Delta }}_1{I}^{\ast }},\\ {\mathrm{{\rm Y}}}_d^{G^{\ast }}& =\frac{d}{G^{\ast }}\frac{\mathrm{\partial }{G}^{\ast }}{\mathrm{\partial }d}=\frac{bd{I}^{\ast }}{{\mathrm{\Delta }}_1(e+I^{\ast })},{\mathrm{{\rm Y}}}_d^{I^{\ast }}=\frac{d}{I^{\ast }}\frac{\mathrm{\partial }{I}^{\ast }}{\mathrm{\partial }d}=-\frac{d(a+b{I}^{\ast })}{{\mathrm{\Delta }}_1(e+I^{\ast })},\end{array}$where ${\mathrm{\Delta }}_1=bc{G}^{\ast }+\frac{de}{(e+I^{\ast }{)}^2}(a+b{I}^{\ast })>0$. By direct calculation, it is easy to obtain $\begin{array}{rl}|{\mathrm{{\rm Y}}}_c^{G^{\ast }}|& =|{\mathrm{{\rm Y}}}_b^{I^{\ast }}|<1,\\ |{\mathrm{{\rm Y}}}_b^{G^{\ast }}|& <|{\mathrm{{\rm Y}}}_d^{G^{\ast }}|<|{\mathrm{{\rm Y}}}_d^{I^{\ast }}|,\\ |{\mathrm{{\rm Y}}}_c^{G^{\ast }}|& =|{\mathrm{{\rm Y}}}_b^{I^{\ast }}|<|{\mathrm{{\rm Y}}}_c^{I^{\ast }}|.\end{array}$The results show that the size of the increase in $c(b)$ will lead to a smaller size of the decrease in $G^{\ast }(I^{\ast })$. $G^{\ast }$ is more sensitive to change d than b, and $I^{\ast }$ is more sensitive to change c than b. In view of the absolute values, we can get the size of changing $I^{\ast }$ is larger than $G^{\ast }$ with changing the same proportion of d as well as c.

6. The strategies and numerical stimulations

The purpose of this paper is to attempt to control the glucose level G in normal range by changing some parameters, and we hope to find some treatments which could be applied clinically according to our model. For diabetic patients, the plasma glucose is higher than healthy persons, then a natural advice is to eat less carbohydrates, i.e. the $G_{in}$ is decreased. Relative to insulin-independent utilization rate a, we are more concerned about other parameters b, c and d, which can affect the state of glucose level. Then we always fix parameter a as a constant. In addition, we are also concerned about the effect of insulin dose on blood glucose. In the following, we always adopt the glucose standards announced by WHO (1999) that the value of normal glucose range is $70<G<110(\mathrm{mg}/\mathrm{dl})$ after fasting and below $140(\mathrm{mg}/\mathrm{dl})$ after postprandial 2 hours [20]. To illustrate the theoretical results, some numerical simulations are displayed in Figures 2 and 3. Some of the values of parameters are taken in the reference [10,15], and the units of the parameters are converted as needed. Four strategies are designed as follows.

Figure 2. The glucose level is affected by insulin dose injected, insulin resistance, insulin secretion and clearance rate. (a) The glucose level changes with the insulin injection when fixing parameters $b=7.5919\times {10}^{-4},c=0.2298,d=1500$. (b) The glucose level changes with the insulin resistance if no insulin injection and c = 0.2298, d = 1500. (c) The glucose level changes with the insulin secretion if no insulin injection and $b=7.5919\times {10}^{-4},d=1500$. (d) The glucose level changes with the insulin clearance rate if no insulin injection and $b=7.5919\times {10}^{-4},c=0.2298$.

Figure 3. Solution trajectories of the model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}G}{\mathrm{d}t}=G_{in}-aG-bGI,}\\ {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=I_{in}+cG-\frac{dI}{e+I},}\end{array}$(1) ) with respect to the subjects data [15]. (a) The glucose level changes with different dose insulin injection for subject Female 3: $a=0.0002,b=7.5919\times {10}^{-4},c=0.2298$. (b) The glucose level changes with different dose insulin for mean of five subjects: $a=0.0309,b=1.6089\times {10}^{-4},c=0.452725$.

Strategy 1: Adjust the dose of insulin injections.

For T1DM, because their damaged β cells secret insulin hardly, the artificial insulin injection will stay with them for the rest of their lives to help decomposing glucose, as well as some T2DM are also recommended to inject insulin in current cure process. So we first consider how to adjust the dose of insulin injection to obtain the normal glucose range. For fixed $G_{in},a,b,c,d,e$, we can obtain $I^{\ast }$ from first equation of (2(2) $\{\begin{array}{l}{G}_{in}-aG-bGI=0,\\ {\displaystyle I_{in}+cG-\frac{dI}{e+I}=0.}\end{array}$(2) ) and substitute it into second equation of (2(2) $\{\begin{array}{l}{G}_{in}-aG-bGI=0,\\ {\displaystyle I_{in}+cG-\frac{dI}{e+I}=0.}\end{array}$(2) ), we have $I_{in}=\frac{d\frac{G_{in}-a{G}^{\ast }}{b{G}^{\ast }}}{e+\frac{G_{in}-a{G}^{\ast }}{b{G}^{\ast }}}-c{G}^{\ast }=\frac{d(G_{in}-a{G}^{\ast })}{eb{G}^{\ast }+(G_{in}-a{G}^{\ast })}-c{G}^{\ast }.$The result indicates that $G^{\ast }$ is able to keep in normal range by adjusting the injection dose of $I_{in}$ (see Figure 2 a). The profile shows $G^{\ast }$ decreases as the insulin dosage $I_{in}$ increases.

For T2DM, their β cells secret insulin, but there are two cases that the glucose cannot be broken down properly. One is the secreted insulin is not enough, another is the insulin resistance occurs unusually through there are enough insulin (even excess). Therefore, under the condition of no insulin injection, we focus on the situation and design three strategies for T2DM.

Strategy 2: Control the insulin resistance.

Suppose there is no insulin injection, i.e. $I_{in}=0$, for fixed $a,c,d,e,G_{in}$, it is easy to see $b=\frac{G_{in}-a{G}^{\ast }}{G^{\ast }{I}^{\ast }}=\frac{G_{in}-a{G}^{\ast }}{G^{\ast }\frac{ce{G}^{\ast }}{d-c{G}^{\ast }}}=\frac{(G_{in}-a{G}^{\ast })(d-c{G}^{\ast })}{ce{G^{\ast }}^2}.$In combination with the relationship of b and $G^{\ast }$ in Section 5, it indicates that improving the ability of insulin resistance b could help decreasing the concentration of $G^{\ast }$(see Figure 2 b). Currently, the main drug used to improve insulin sensitivity is metformin and analogue.

Strategy 3: Stimulate insulin secretion.

Suppose there is no insulin injection, i.e. $I_{in}=0$, for fixed $a,b,d,e,G_{in}$, we have $c=\frac{\frac{d+I^{\ast }}{e+I^{\ast }}}{G^{\ast }}=\frac{d\frac{G_{in}-a{G}^{\ast }}{b{G}^{\ast }}}{G^{\ast }(e+\frac{G_{in}-a{G}^{\ast }}{b{G}^{\ast }})}=\frac{d(G_{in}-a{G}^{\ast })}{G^{\ast }(eb{G}^{\ast }+G_{in}-a{G}^{\ast })}.$The structure of the formula shows the insulin secretion rate c could affect the glucose $G^{\ast }$. If the rate c can be adjusted, then the level of $G^{\ast }$ can be controlled in normal value (see Figure 2 c). Glimepiride is used to stimulate insulin secretion widely, as well as gliclazide and so on.

Strategy 4: Adjust the insulin clearance rate.

From the second equation of (2(2) $\{\begin{array}{l}{G}_{in}-aG-bGI=0,\\ {\displaystyle I_{in}+cG-\frac{dI}{e+I}=0.}\end{array}$(2) ), we have $(I_{in}+c{G}^{\ast })(e+I^{\ast })=d{I}^{\ast }$. Still consider there is no insulin injection, i.e. $I_{in}=0$, and substitute it into the first equation of (2(2) $\{\begin{array}{l}{G}_{in}-aG-bGI=0,\\ {\displaystyle I_{in}+cG-\frac{dI}{e+I}=0.}\end{array}$(2) ), it is easy to get $d=\frac{(I_{in}+c{G}^{\ast })(e+\frac{G_{in}-a{G}^{\ast }}{b{G}^{\ast }})}{\frac{G_{in}-a{G}^{\ast }}{b{G}^{\ast }}}=\frac{c{G}^{\ast }(be{G}^{\ast }+G_{in}-a{G}^{\ast })}{G_{in}-a{G}^{\ast }}.$

The result indicates that $G^{\ast }$ can be controlled in normal value by adjusting the maximum insulin clearance rate d (see Figure 2 d).

To illustrate the effects of dosage insulin injected on plasma glucose, we take the data from subject Female 3 and mean of five subjects for example [15]. We set $d=1500,e=2300,G_{in}=4.5$. In Figure 3(a), take the data of subject Female 3 $a=0.0002,b=7.5919\times {10}^{-4},c=0.2298$, and give $I_{in}=0,5,10,20$, then the corresponding $G^{\ast }$ is 128.052, 117.526, 107.912, 91.2907. Obviously, if there is no insulin injected, $G^{\ast }=128.052$ is above the normal glucose. Once given insulin, the level of glucose will be controlled in proper range. In Figure 3(b), we take $a=0.0309,b=1.6089\times {10}^{-4},c=0.452725$, and still give $I_{in}=0,5,10,20$ respectively, then the corresponding $G^{\ast }$ is $116.708,112.272,108.174,100.834$. To sum up, the insulin plays an important part in decreasing the level of glucose. Clinically, measuring diabetics patient's insulin index value and determining the dose of extracorporeal insulin will contribute to adjust glucose more precisely.

Remark: For the patients with T1DM, they only have a quite small portion of intact β cells, which lead to very little insulin secreted. Then it is reasonable to assume the insulin secretion rate c is close to zero. Based on this case, to keep glucose G in the normal range, insulin injection has become the preferred therapy for T1DM so far clinically. For T2DM, the insulin secreted by β cells is not enough to help decomposing glucose, so the secretion rate $c\in (0,1]$ [10,21]. Under the circumstances, the early T2DM patients can take some medicines such as melbine, gliclazide and glimepiride and so on, which is able to help secreting insulin or increasing insulin sensitivity. Once those medicines are out of action, they need to inject insulin for controlling glucose in the desirable range.

7. Summary and discussion

For diabetics patients, the main aim is to adjust the glucose level in normal range without changing the life style of patients. As well known, the current therapies are mainly insulin injections and oral hypoglycemic drugs, such as metformin, gliclazide and so on. The treatment mechanisms of these medicines are different. With the help of artificial insulin, we can control the glucose level easily, but hypoglycemia and hyperglycemia appear if the dose of insulin is not right. Therefore, it is very important that how to regulate insulin better in the treatment of diabetes. Undoubtedly, patients with type 1 diabetes need to inject insulin daily because of deficient insulin production absolutely. The insulin pump is the most advanced way of insulin administration for T1DM currently. Patients with type 2 diabetes are called non-insulin-dependent because the pancreas does not produce enough insulin or the body cannot use the secreted insulin effectively. For early T2DM, some of them are always advised to inject insulin in current clinical practice. Is the insulin injection unique and best way? Actually, lots of experiments show that the level of plasma glucose could be affected by individual diabetics patient factors, such as weight, age, food-intake, exercise and sensitivity to insulin. In addition, for individual patient, the time taken to absorb one-half of an injected dose of insulin also may change by 25–50% [22]. For them, if there are some treatments that do not use insulin injection, then the financial burden and physiological pain are both alleviated. In [10], authors considered $G_{in}$ and $I_{in}$ were positive periodic functions and mimicked pancreatic insulin secretion by exogenous insulin infusions. Their results illustrated the preprandial insulin could better control the postprandial glucose level. But they did not mention other factors affecting glucose control. Nevertheless, who are suitable for insulin injections and who are suitable for oral medicines, there is no uniform standard. Therefore, it is especially important to design personalized treatments for different individuals. With proper insulin treatments, patients can maintain plasma glucose within a normal range but hypoglycemia and hyperglycemia could not happen. In the paper, we attempt to design optimal therapies which could help the patients to reduce pain and achieve the best hypoglycemic effect. A glucose–insulin regulatory model (1(1) $\{\begin{array}{l}{\displaystyle \frac{\mathrm{d}G}{\mathrm{d}t}=G_{in}-aG-bGI,}\\ {\displaystyle \frac{\mathrm{d}I}{\mathrm{d}t}=I_{in}+cG-\frac{dI}{e+I},}\end{array}$(1) ) with Michaelis–Menten function as insulin degradation rate is studied analytically and numerically. We have shown that it is generally possible to individually alter these parameters b, c and d to stabilize average glucose levels except injecting insulin $I_{in}$ for T2DM. Therefore, they can avoid the pain of subcutaneous insulin injections daily before their β cells loss functions completely. According to different situations, the doctors could choose one of above four strategies or a combination of treatments discussed in Section 6 for serious diabetic patients. We also point out that the liver provides glucose when the hypoglycemia happens. This is a more complex dynamic which is worth studying in our future work.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work is supported by the Ph. D. Fellowship of Xi'an Polytechnic University (No. 107020344), The Natural Science Basic Research Plan in Shaanxi Province of China (No. 2021JM-444), CSC (No.201908610049).