Dynamic analysis of a Filippov blood glucose insulin model

1.   Introduction

Diabetes mellitus is a chronic metabolic disease characterized by hyperglycemia and hypoglycemia, which is mainly divided into four types; these include type 1 diabetes, type 2 diabetes, gestational diabetes, and other types of diabetes; it can cause severe complications ^[1]. Due to the damage of $\beta$ cells in the patient's pancreas, they secrete almost no insulin or can only secrete a small amount of insulin, and the patient can only rely on subcutaneous injections of exogenous insulin to regulate their blood glucose concentration. However, due to the dysfunction of the blood glucose insulin regulatory system, so-called insulin resistance is caused, which makes insulin unable to be effectively utilized. Therefore, the cells need to synthesize and secrete more insulin to eliminate insulin resistance ^[2]. The typical clinical manifestations of diabetes are polydipsia, polyuria, and weight loss. A long course of this disease can lead to chronic progressive lesions of many tissues and organs in the body, such as eyes, kidneys, nerves, and the heart. In the past decades, the prevalence and incidence rate of diabetes has been rising. Because diabetes needs lifelong treatment, it causes high disease burden and economic burden to patients and families, and it has become a major public health problem in the world.

The main treatment of diabetes is to control the blood sugar concentration within a certain range through the input of exogenous insulin and its analogues; to reduce and delay the occurrence of complications. Insulin therapy is an effective method for diabetes patients to control hyperglycemia, as it can reduce blood sugar and promote the synthesis of glycogen, fat, and protein. Both insulin pumps and artificial pancreases can be used for the treatment of diabetes. Insulin pumps are the most common medical device for administering insulin and its analogues. It can enable patients to decide the injection dose and injection time according to their needs, but it is easy to cause various complications due to poor blood glucose control ^{[3,4,5,6,7,8]}. An artificial pancreas is a device that is able to analyze information from sensors, such as continuous glucose monitoring, to deliver the correct amount of insulin by subcutaneous injection via a pump; it can exert the endocrine function of the pancreas in diabetes, and automatically control the glucose concentration in the body within the normal range ^[9,10,11]. Whether it is an insulin pump or an artificial pancreas, the estimation and distribution of insulin injection dose should follow different models and rules according to the type of diabetes, blood sugar level, and individual conditions of the patient, and personalized dose adjustment should be made according to the blood sugar detection level during use ^{[12,13,14,15]}. At present, the most comprehensive insulin treatment measures have not been sourced. Therefore, there is a lot of research space for insulin treatment measures.

In order to explore reasonable and effective treatment methods for controlling blood glucose concentration, many researchers have studied the mechanism of blood glucose insulin interaction by establishing relevant mathematical models. Bolie et al. first proposed a blood glucose insulin interaction model in 1961 to simulate the sub-diurnal oscillation process of insulin secretion and the regulatory mechanism of the metabolic system in the human body ^[16]. Considering the function of $\beta$-cell in maintaining blood insulin concentration, Abdul et al. proposed a new time-delay mathematical model for the glucose-insulin endocrine metabolism regulation feedback system ^[17,18]. Angelova et al. considered the glucose insulin interaction process in the human body and established a two-dimensional delayed differential equation system ^[19]. Muhammad et al. proposed an impulsive differential equation model to study plasma glucose control for diabetic patients with impulsive insulin injections ^[20]. Iqra et al. introduced an improved Hovorka model to further explain glucose-insulin kinetics ^[21]. Li et al. established an interaction model to study the relationship between changes in blood glucose concentration and insulin content by introducing two time delays ^[22]. Ling et al. established an ordinary differential equation model with three physiological delays by taking into account the transport delay of insulin from the secretion site to the plasma ^[23]. Ma et al. proposed a glucose-insulin model with Michaelis-Menten function as the insulin degradation rate to mimic the pathogenic mechanism of diabetes ^[24]. Rao et al. investigated the dynamics of a glucose-insulin regulatory system model that considers discrete time delays ^[25]. Shi et al. proposed a glucose-insulin regulatory system with the perturbation of white noise and the mass of $\beta$-cell as variables ^[26,27]. The above studies provide guidance for clinical doctors to control blood sugar concentration within the ideal range.

Because diabetes is incurable, diabetes treatment is a long-term and uninterrupted process. Based on the characteristics of insulin injection behavior, the use of pulse differential equations to describe the dynamic relationship between glucose and insulin has been widely adopted by many scholars. However, many pulse systems cannot describe the process of blood glucose growth and insulin treatment ^[9,28,29,30]. Filippov system is a discontinuous dynamical system composed of two or more smooth vector fields, where different smooth vector fields are divided by discontinuous boundaries. In fact, this system has many important applications in biological, neuronal, and mechanical fields ^[31,32,33]. Therefore, establishing a Filippov model with a threshold strategy ^[34,35] is of great practical significance.

This paper is arranged as follows: In Sections 2 and 3, a Filippov model is established to analyze the dynamic properties of the two subsystems. In Section 4, we discuss the stability of the sliding line region and equilibrium point of the system. In Section 5, we validate the theoretical analysis results through numerical simulations. In Section 6, we briefly discuss and summarize the main results.

2.   Model formulation

At present, researchers are mainly dedicated to the development of artificial pancreases, which consists of several glucose concentration monitoring sensors, an insulin injection device and a control algorithm. Artificial pancreases can continuously monitor the glucose concentration in the blood and provide the corresponding amount of insulin. The glucose concentration monitoring system obtains more than 200 blood glucose measurements every day, which can obtain the trend of blood glucose concentration changes and make timely adjustments. Therefore, open-loop control can be transformed into closed-loop control through a glucose monitoring system with single or multiple glucose sensors. When the blood glucose concentration feedback from the blood glucose monitoring system reaches a specific threshold, insulin injection can be performed. Injecting an appropriate amount of insulin into the body can greatly reduce the occurrence of hypoglycemia and hyperglycemia. Considering the working principle of the artificial pancreas, Huang et al. constructed a semi-continuous dynamic system with state feedback insulin pulse injection ^[9]:

where $G(t)$ and $I(t)$ represents the glucose concentration and insulin concentration at time $t\geq 0$, respectively. $G_{in}$ represent the constant glucose exogenous infusion, $\sigma_2 G$ and $a(c+ \frac{mI}{n+I})G$ stand for the insulin-independent and insulin-dependent glucose consumption, respectively, and $b$ is the glucose production rate. $\frac{\sigma_1 G^2}{\alpha_1^2 + G^2}$ is the insulin secretion stimulated by elevated glucose concentration, $d_i I$ is the constant degradation of insulin, and $\sigma$ is the dose in each injection. Here, $\sigma$, $\sigma_1$, $\sigma_2$, $\alpha_1$, $m$, $n$, $a$, $b$, $c$, and $d_i$ are positive constant parameters. $M_0 = \left\{(G, \; I)|G = E_G \text{ and } I > I_G\right\}$, $E_G$ is a preset threshold, and $I_G$ can be obtained through model parameter calculation. This model analyzes the controllability of blood glucose levels by closely monitoring blood glucose levels, i.e., insulin injection when blood glucose levels reach or exceed preset adjustable thresholds.

Based on the above model, this paper considers how to maintain a blood glucose concentration level through external intervention. It is possible to set a blood glucose level for external intervention through methods such as blood glucose monitoring systems. When the blood glucose concentration drops below a specific control threshold $E_G$, it needs to suspend insulin injection or reduce the insulin injection dose; when the blood glucose concentration reaches the control threshold $E_G$, insulin injection is needed to reduce blood glucose.

Due to the sub-diurnal oscillation process of insulin secretion and the regulatory mechanism of the human metabolic system, we chose to ignore the physiological delays in insulin secretion, transportation, and liver glycogen breakdown. We set insulin secretion as a constant when glucose concentration increases. When the blood glucose concentration is below the control threshold $E_G$, no treatment measure is required. Here, the mechanism of blood glucose insulin interaction can be described as:

where $\sigma_1$ is the insulin secretion constant when glucose concentration increases.

When the blood glucose concentration reaches or exceeds the control threshold $G = E_G$, a therapeutic measure is required, such as injecting insulin, to promote a decrease in blood glucose concentration. Now, the consumption of insulin dependent sugars will increase, so the model is transformed into

where $c_1 G$ and $\sigma$ represent glucose consumption and insulin injection dependent on insulin after treatment measure is taken, respectively.

Let $z = (G, \; I)^T, \; H(z) = G-E_G,$ and

The following Filippov model can be obtained:

where

Systems (2.1) and (2.2) are subsystems of system (2.3). Also, $\Sigma = \{z\in\mathrm{R}^2_+:H(z) = 0\}$ is denoted as the boundary of separating the two regions $\mathbb{S}_1$ and $\mathbb{S}_2$, and $H$ is a smooth scalar function with nonzero gradients on $\Sigma$ and $H_z = (1, 0)$.

3.   The dynamic properties of subsystem

3.1.   The dynamic properties of subsystem (2.1)

System (2.1) can be abbreviated as

where $\mathbb{G} = G_{in}+b, \sigma_3 = \sigma_2+ac$.

(1) When $\sigma_1 = 0$, it is type 1 diabetes, and the $\beta$ cells in the patient's own pancreas secrete little or no insulin; thus, the patient can only rely on the subcutaneous injection of exogenous insulin to regulate their blood sugar concentration. System (3.1) has a unique equilibrium point $E_0 = (G^0, 0)$, where $G^0 = \frac{\mathbb{G}}{\sigma_3}$.

Calculating the Jacobian matrix of system (3.1) at $E_0$, we obtain that

Its characteristic polynomial is $(\lambda+\sigma_3)(\lambda+d_i) = 0$, where two characteristic roots are $\lambda_1 = -\sigma_3 < 0$, $\lambda_2 = -d_i < 0$. Therefore, $E_0$ is a locally asymptotically stable node.

(2) When $\sigma_1 > 0$, it is type 2 diabetes, and the $\beta$ cells in the patient's body can secrete a small amount of insulin. System (3.1) has a unique positive equilibrium point $E_1 = (G_1, I_1)$, where

Theorem 3.1. System (3.1) has a unique positive equilibrium point $E_1$ when $\sigma_1 > 0$, and it is globally asymptotically stable.

Proof. Let the right end of the second equation of system (3.1) be $Q_1 = 0$, and we get by calculation that $G = d_1I$. Inserting $G = d_1I$ into equation $P_1 = 0$, we have that

Obviously, $-n\mathbb{G} < 0$. Notice the discriminant of quadratic polynomials:

Therefore, the above univariate quadratic equation has a unique positive root, and let $I_1 = \frac{(\mathbb{G}-n\sigma_3 d_1+\sqrt{\Delta_1})\sigma_1}{2d_1(\sigma_3+am)}$. Subsequently, we get that the system (3.1) has a unique positive equilibrium point $E_1$.

Calculating the Jacobian matrix of the system (3.1) at $E_1$, we have that

Its characteristic polynomial is $\lambda^2+a_1\lambda+a_2 = 0$, where

Therefore, $E_1$ is a locally asymptotically stable node or focus point.

In addition, define the Dulac function $B(G, \; I) = \frac{1}{GI}$ on a simply connected region $\mathrm{R}^2_+$, then we have that

According to Dulac's theorem, system (3.1) has no limit cycle within $\mathrm{R}^2_+$. Therefore, $E_1$ is globally asymptotically stable on $\mathrm{R}^2_+$. This completes the proof.

3.2.   The dynamic properties of subsystem (2.2)

System (2.2) can be rewritten as

where $\sigma_4 = \sigma_3+ac_1$.

We will prove that system (3.2) has a unique positive equilibrium point $E_2 = (G_2, \; I_2)$, where

and $\Delta_2 = (\sigma am+\sigma_1\mathbb{G}+\sigma\sigma_4-\sigma_4nd)^2+4nd_i(\sigma_4+am)(\sigma\sigma_4+\sigma_1\mathbb{G})$.

Theorem 3.2. The system (3.2) has a unique positive equilibrium point $E_2$ and it is globally asymptotically stable.

Proof. Let $Q_2 = 0$, and we get that $G = \frac{d_iI-\sigma}{\sigma_1}$. Substituting $G$ into $P_2 = 0$ yields

Its discriminant is

Further, we get the unique positive root, denoted as

Thus, we have that

If $G_2 > 0$, we have

By calculations, we get that $\sqrt{\Delta_2} > 2\sigma(\sigma_4+am)-(\sigma am+\sigma_1\mathbb{G}+\sigma\sigma_4-\sigma_4nd)$. Squaring both sides yields that $dn+\sigma > 0$ must hold, thus $G_2 > 0$.

Calculating the Jacobian matrix of system (3.2) at $E_2$, we have

Its characteristic polynomial is $\lambda^2+a_3\lambda+a_4 = 0$, where

Therefore, $H_1 = a_3 > 0, \; H_2 = a_3a_4 > 0$. According to Hurwitz criterion, we have that $E_2$ is a locally asymptotically stable node or focus point.

Obviously, $\mathrm{R}^2_+$ is a simply connected region. Defining the Dulac function $B(G, \; I) = \frac{1}{GI}$, we have

According to Dulac's theorem, system (3.2) has no limit cycle within $\mathrm{R}^2_+$. Therefore, $E_2$ is globally asymptotically stable on $\mathrm{R}^2_+$. This completes the proof.

4.   Sliding region and equilibrium point of Filippov system

This section discusses the sliding region and equilibrium point of system (2.3). To begin, we give the definitions of the sliding region and equilibrium point.

Consider a class of Planar Filippov system ^[36]

where $\mathbb{S}_1 = \{z\in R^+_2:\; H(z) < 0\}, \; \mathbb{S}_2 = \{z\in R^+_2:\; H(z) > 0\}$, and $\; R^+_2 = \{z = (x, y:x\geq0, y\geq0)\}.$ In addition, define $\Sigma = \{z\in\mathrm{R}^2_+:H(z) = 0\}$ as the switching surface of separating the two regions $\mathbb{S}_1$ and $\mathbb{S}_2$, and $H$ is a smooth scalar function with nonzero gradients on $\Sigma$ and $H_z = (1, 0)$.

Definition 4.1. (^[36,37]) $\Sigma_s = \{z\in\Sigma|\sigma_z\leq 0\}\subseteq\Sigma$ is called a sliding region of system (2.3) if $\sigma_z = \langle H_z, \; F_{\mathbb{S}_1}(z)\rangle \cdot \langle H_z, \; F_{\mathbb{S}_2}(z)\rangle = F_{\mathbb{S}_1}(z)\; H_z \cdot F_{\mathbb{S}_2}(z)\; H_z$, where $< , >$ denotes the standed scalar product of $H(z)$ and $F_{\mathbb{S}_i}, i = 1, 2$.

$\Sigma_s$ is stable if $\langle H_z, \; F_{\mathbb{S}_1}(z)\rangle > 0$ and $\langle H_z, \; F_{\mathbb{S}_2}(z)\rangle < 0$, whereas it is unstable if $\langle H_z, \; F_{\mathbb{S}_1}(z)\rangle < 0$ and $\langle H_z, \; F_{\mathbb{S}_2}(z)\rangle > 0$. In particular, $z\in\Sigma_s$ is called a singular sliding point if $\sigma_z = 0$.

Definition 4.2. (^[36,37,38]) $z$ is called a regular equilibrium point of system (2.3) (denoted as $E^R$) if $z\in \mathbb{S}_1$, $F_{\mathbb{S}_1}(z) = 0$ or $z\in \mathbb{S}_2$, $F_{\mathbb{S}_2}(z) = 0$.

$z$ is called a virtual equilibrium point of system (2.3) (denoted as $E^V$) if $z\in \mathbb{S}_1$, $F_{\mathbb{S}_2}(z) = 0$ or $z\in \mathbb{S}_2$, $F_{\mathbb{S}_1}(z) = 0$.

$z$ is called a boundary equilibrium point of system (2.3) (denoted as $E^B$) if $F_{\mathbb{S}_1}(z) = 0, H(z) = 0$ or $F_{\mathbb{S}_2}(z) = 0, H(z) = 0$.

$z$ is called a tangent point of system (2.3) (denoted as $E^T$) if $z\in\Sigma_s, \; F_S(z)\neq 0$ and $< F_{\mathbb{S}_1}\cdot H(z) > < F_{\mathbb{S}_2}\cdot H(z)\geq0$.

Definition 4.3. (^[36,37]) The local trajectory on $\Sigma_s$ is defined by the convex combination of the Filippov system. Considering system

where $z\in\Sigma_s, \; \Sigma_s(z)$ is called the sliding system of system (2.3), $z$ is a pseudo equilibrium point (denoted as $E^p$) if $Z_s(z) = 0$.

Obviously, we know that $G_2 < G_1$, thus the authenticity of the positive equilibrium point of system (2.3) can be addressed in three different scenarios:

(1) If $\frac{\mathbb{G}-n\sigma_3d_1+\sqrt{\Delta_1}}{2(\sigma_3+am)} < E_G < \frac{\sigma am+\sigma_1\mathbb{G}+\sigma\sigma_4-\sigma_4nd+\sqrt{\Delta_2}}{2(\sigma_4+am)}-\frac{\sigma}{\sigma_1}$, then $E_1$ is a virtual equilibrium point (denoted as $E_1^V$), and $E_2$ is a regular equilibrium point (denoted as $E_2^R$).

(2) If $\frac{\sigma am+\sigma_1\mathbb{G}+\sigma\sigma_4-\sigma_4nd+\sqrt{\Delta_2}}{2(\sigma_4+am)}-\frac{\sigma}{\sigma_1} < E_G < \frac{\mathbb{G}-n\sigma_3 d_1+\sqrt{\Delta_1}}{2(\sigma_3+am)}$, then $E_1, \; E_2$ are both virtual equilibriums.

(3) If $E_G > \frac{\mathbb{G}-n\sigma_3 d_1+\sqrt{\Delta_1}}{2(\sigma_3+am)}$, then $E_1$ is a regular equilibrium point (denoted as $E_1^R$), and $E_2$ is a virtual equilibrium point (denoted as $E_2^V$).

Further, we discuss the existence of a sliding region. According to the Definition 4.1, we get the sliding region $\Sigma_s = \{z\in \Sigma|\sigma(z)\leq0\}$.

Theorem 4.1. If $E_G < \frac{G_{in}+b}{\sigma_2+ac}$, system (2.3) has a sliding region.

Proof. By calculation, we have

Solving the inequality $\sigma(z)\leq 0$ yields

where

The sliding region of system (2.3) exists, i.e.,

To ensure that the existence of the sliding region is always true, we assume that $G_{in}+b-(\sigma_2+ac)E_G > 0$; and we get that $E_G < \frac{G_{in}+b}{\sigma_2+ac}$. This completes the proof.

Furthermore, according to the Definition 4.2, we have that the systems (2.1) and (2.2) satisfy $G = E_G$, and there exists a boundary equilibrium point for system (2.3), which is

System (2.3) has a tangent point which satisfies $G = E_G$ in the sliding region $\Sigma_s$. By calculations, we have that

Finally, we now discuss the existence and stability of pseudo equilibrium points.

Theorem 4.2. If $G_2 < E_G < min\{G_1, \; \frac{G_{in}+b}{\sigma_2+ac}\}$, there exists a pseudo equilibrium point and it is locally asymptotically stable.

Proof. According to the Definition 4.3, we have

where

When $Z_s(z) = 0$, i.e., $AI^2+BI+C = 0$. It is assumed that $F(I) = AI^2+BI+C = 0$, due to $C < 0$. By calculations, the discriminant $\Delta_3 > 0$ is always true, and the unique positive solution can be obtained by

where

If there exists a pseudo equilibrium point, then it needs to meet the following requirements: $I \in \Sigma_s$ i.e., $\{I\; |\; I_1\leq I\leq I_2\}$. In this case, we have a pseudo equilibrium point, denoted as $E_p(E_G, I_p)$.

Due to $F'(I_p) = -\sqrt{\Delta_3} < 0$, it can be concluded that the pseudo equilibrium point $E_p$ is locally asymptotically stable. This completes the proof.

Based on the above discussion, the dynamic behavior of system (2.3) can be obtained that:

(1) If $E_G < E_2 < E_1$, the regular equilibrium point $E_2^R$ of system (2.3) is globally asymptotically stable;

(2) If $E_G = E_2$, system (2.3) has boundary equilibrium points $E_2^B$ and virtual equilibrium points $E_1^V, E_2^V$;

(3) If $E_2 < E_G < \min\{E_1, \frac{G_{in}+b}{\sigma_2+ac}\}$, the pseudo equilibrium point $E_p$ of system (2.3) is locally asymptotically stable;

(4) If $E_G = E_1$, system (2.3) has boundary equilibrium points $E_1^B$ and virtual equilibrium points $E_1^V, E_2^V$;

(5) If $E_2 < E_1 < E_G$, the regular equilibrium point $E_1^R$ of system (2.3) is globally asymptotically stable.

5.   Numerical simulations

This section shows numerical simulations on different types of equilibrium points for system (2.3) to confirm the above theoretical results. For sake of convenience, we have used the same approach as ^[9,14,39] to convert units from mass to concentration.

5.1.   Numerical example 1

The parameter values in our simulations are the same as in ^[10,19]. We select parameters $G_{in} = 216 \; \rm mg/min, \; m = 900 \; \rm mg/min, \; n = 80 \; \rm mg$, $a = 0.00003 \; \rm mg^{-1}, \; b = 100 \; \rm mg/min$, ${\rm{c}} = 40 \; \rm mg/min, \; c_1 = 30 \; \rm mg/min$, $\sigma = 20 \; \rm mg, \; \sigma_1 = 0.002 \; \rm min^{-1}$, $\sigma_2 = 0.000005 \; \rm min^{-1}, \; d_i = 0.08 \; \rm min^{-1}.$ By calculations, we get that the positive equilibrium points of the two subsystems are $E_1 = (G_1, I_1) = (136.8648, 34.216), \; E_2 = (G_2, I_2) = (122.8844, 55.721).$ There are three types of equilibrium points that can be determined for $\rm Filippov$ system (2.3):

(1) If $E_G < G_2$, system (2.3) has a regular equilibrium point $E_2$;

(2) If $G_2 < E_G < G_1$, system (2.3) has a pseudo equilibrium point $E_p$;

(3) If $G_1 < E_G$, system (2.3) has a regular equilibrium point $E_1$.

The above analytical results will be validated using numerical simulations. Set the control thresholds for human blood glucose concentration to be $E_G = 110 \; \rm mg/dl$, $\; E_G = 130 \; \rm mg/dl\;$ and $\; E_G = 150 \; \rm mg/dl$, respectively. We have that:

(1) When $E_G = 110 \; \rm mg/dl$, $E_2$ is a regular equilibrium point; and $E_1$ is a virtual equilibrium point. According to Theorem 3.2, the equilibrium point $E_2$ is globally asymptotically stable, and the regular equilibrium point $E_2^R$ is a globally asymptotically stable node (see Figure 1(a)).

(2) When $E_G = 130 \; \rm mg/dl$, $E_1, E_2$ are both virtual equilibrium points; and $E_p$ is a pseudo equilibrium point. According to Theorem 4.2, the pseudo equilibrium point $E_p$ is stable. Therefore, all trajectories ultimately tend toward the pseudo equilibrium point $E_p$ (see Figure 1(b)).

(3) When $E_G = 150 \; \rm mg/dl$, $E_1$ is a regular equilibrium point; and $E_2$ is a virtual equilibrium point. According to Theorem 3.1, the equilibrium point $E_1$ is globally asymptotically stable. At this point, the regular equilibrium point $E_1^R$ is a globally asymptotically stable node (see Figure 1(c)).

Through analysis, it was found that: (1) When $E_G < G_2$, the blood glucose concentration in the human body tends to be normal, and no treatment measures need to be taken (i.e., stopping insulin injection or reducing insulin injection dose); (2) when $G_2 < E_G < G_1$, the pseudo equilibrium point always exists and is stable, and the plasma glucose concentration in the human body is higher than the control threshold. Treatment measures need to be taken, and an appropriate amount of insulin injection can be performed to control blood sugar at a certain level; (3) when $G_1 < E_G$, the blood glucose concentration will not be constrained by threshold control strategies and insulin therapy is needed. This indicates that when the human blood glucose concentration control threshold is set reasonably, the blood glucose concentration can be stably controlled within the normal range.

5.2.   Numerical example 2

In this example, we replace the $\sigma_1$ (i.e., insulin secretion when glucose concentration increases) in system (2.3) by a functional reflection function $f = \frac{\sigma_1G^2}{\alpha_1^2+G^2}$ to better fit the specific geometric characteristics exhibited in laboratory data fitting.

Here, system (2.1) is converted as

and system (2.2) is converted as

thus, $\rm Filippov$ system (2.3) is converted as

where

We select parameters $G_{in} = 216 \; \rm mg/min, \; m = 900 \; \rm mg/min, \; n = 80 \; \rm mg$, $a = 0.00003 \; \rm mg^{-1}, \; b = 100 \; \rm mg/min$, ${\rm{c}} = 40 \; \rm mg/min, \; c_1 = 30 \; \rm mg/min$, $\sigma_2 = 0.000005 \; \rm min^{-1}, \; \sigma_1 = 6.27 \; \rm mU/min$, $\sigma = 130 \; \rm mU, \; d_i = 0.08 \; \rm min^{-1}$, $\alpha_1 = 105 \; \rm mg.$ By calculations, we have that the positive equilibrium points of systems (5.1) and (5.2) are $E_3 = (G_3, I_3) = (216.9387, 7.837), \; E_4 = (G_4, I_4) = (132.6857, 32.837)$. When $E_G$ increases, the equilibrium point type of $\rm Filippov$ system (5.3) exhibits the following pattern :

(1) The virtual equilibrium point $E_3^V$ and the regular equilibrium point $E_4^R$ coexist;

(2) The boundary equilibrium point $E_4^B$ exists;

(3) The virtual equilibrium point $E_3^V, \; E_4^V$ and sliding region $\Sigma_s$ coexist;

(4) The boundary equilibrium point $E_3^B$ exists;

(5) The regular equilibrium point $E_3^R$ and the virtual equilibrium point $E_4^V$ coexist.

It is clear that bifurcations can occur when the parameter $E_G$ varies. The above analytical results will be validated using numerical simulations.

A boundary node bifurcation for system (5.3) may occur when $E_G$ passes through a critical value. Set the control thresholds for human blood glucose concentration to be $E_G = 120 \; \rm mg/dl$, $\; E_G = 180 \; \rm mg/dl\;$ and $\; E_G = 240 \; \rm mg/dl$, respectively. We have that:

(1) When $0 < E_G = 120 \; \rm mg/dl\; < G_4$, then $E_3$ becomes a virtual equilibrium point and $E_4$ is a regular equilibrium point. The equilibrium point $E_4$ is globally asymptotically stable, the regular equilibrium point $E_4^R$ is a globally asymptotically stable node, and $E_3^V$ and $E_4^R$ coexist (see Figure 2(a)). In this case, the blood glucose concentration in the human body tends to be normal; injection of insulin would control the blood glucose concentration at a low level leading to hypoglycaemia, so no treatment measures should to be taken.

(2) When $G_4 < E_G = 180 \; \rm mg/dl\; < G_3$, the virtual equilibrium points $E_3^V, \; E_4^V$ and sliding region $\Sigma_s$ coexist, and a pseudo equilibrium point $E_p^1$ appears. The pseudo equilibrium point $E_p^1$ is stable. Therefore, all trajectories ultimately tend toward the pseudo equilibrium point $E_p^1$ (see Figure 2(b)). This means that the blood glucose concentration stabilizes at the pseudo equilibrium point, which indicates that the blood glucose concentration could be successfully maintained at a desired level within the normal range. Besides, this bifurcation shows how a stable node becomes a stable pseudo equilibrium point.

(3) When $G_3 < E_G = 240 \; \rm mg/dl$, $E_4$ becomes a virtual equilibrium point and $E_3$ is a regular equilibrium point. The equilibrium point $E_3$ is globally asymptotically stable. At this point, the regular equilibrium point $E_3^R$ is a globally asymptotically stable node, and $E_3^R$ and $E_4^V$ coexist (see Figure 2(c)). This implies that the blood glucose concentration stabilizes at a high level, which leads to hyperglycaemia and glucotoxicity. The blood glucose concentration would not be constrained by threshold control strategies and insulin therapy is needed. It is manifested that this bifurcation indicates how a stable pseudo equilibrium point becomes a stable node.

From the above numerical simulation analysis, we find it is possible to make the treatment more reasonable and effective by selecting an appropriate control threshold for human blood glucose concentration in insulin therapy.

6.   Conclusions

This paper proposes a Filippov blood glucose insulin model with threshold control strategy based on the sub diurnal oscillation process of insulin secretion and the regulatory mechanism of the human metabolic system, ignoring the effects of physiological delays such as insulin secretion, transportation, and liver glycogen decomposition. We assume that patients undergo insulin therapy when their blood glucose concentration reaches the control threshold $E_G$. If a patient's blood glucose concentration is below the control threshold $E_G$, the insulin injection dose should be reduced or no treatment measures should be taken.

Using the convex combination definition and related theories of Filippov system, we prove the global stability of its two subsystems. The existence and conditions of the sliding line domain of Filippov system have been given. Different types of equilibrium states of the system have been determined, and the existence and stability of pseudo equilibrium points have been obtained, i.e., when the positive equilibrium points of subsystems (2.1) and (2.2) are both false equilibrium points, selecting the value within the sliding region as the threshold is most appropriate, and there are pseudo equilibrium points in the sliding region, which are stable focal points or nodes.

Through numerical simulations, we have validated the theoretical results of this paper and find an appropriate threshold range, and demonstrate that suitable blood glucose concentration ranges can be found during blood glucose insulin treatments in practical applications, allowing for more scientific and reasonable insulin injections. These results indicate that by selecting an appropriate control threshold and flexibly applying threshold strategies to control the level of blood glucose concentration increase, the human blood glucose concentration can be maintained within an ideal range, thereby achieving a more economical and effective treatment level.

Note that the factors such as insulin sensitivity to glucose tolerance, physiological delay, white noise interference et al. has a certain impact on blood glucose insulin therapy, thus it will be interesting to know how this might affect the dynamics if we included it in our proposed models. Consequently, we plan to address these topics with the aim of improving optimal strategies for the treatment of diabetes in the future research.

Author contributions

Qiongru Wu, Ling Yu, Xuezhi Li, and Wei Li: Establishment of models, theoretical research, numerical simulations; Qiongru Wu: Design and interpretation of the model, theoretical analysis, drafted the manuscript, revision; Wei Li: Design and theoretical analysis of the model, provided critical feedback to the manuscript, revision; Ling Yu: Data collection, numerical simulation, provided suggestions for the manuscript; Xuezhi Li: Revision, translation, proofreading of the manuscript. All authors have read and approved the final version of the manuscript for publication.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this article.

Conflict of interest

The authors declare no conflict of interest.