Glucose-insulin mathematical model for the combined effect of medications and life style of Type 2 diabetic patients

The goal of this paper is to propose a new mathematical model for the combined effect of different treatments and lifestyles on the glucose-insulin dynamics of Type 2 diabetes (T2D) patients. The model gives the possibility to take into consideration physical activity, stress, meals, and medications while evaluating or designing treatment plans for T2D patients. The model is proposed by combining and modifying some of the available models in the literature. Simulations were performed for the modifications to show how the model confirm with literature on T2D patients. Additionally, a discussion is provided to demonstrate the ability of the model to be used in the assessment of treatment plans and in the design for robust insulin dose guidance algorithms. An open source code for the model is additionally provided.


Introduction
One of the greatest health challenges which faces humanity in the 21st century is the emergence of type 2 diabetes (T2D) as a global pandemic. More than 415 million were reported to suffer from T2D in 2015 and the number is expected to reach 642 million by 2040 [5]. Moreover, the global expenses related to T2D are estimated to be 850 million USD in 2017 and they are expected to increase [21]. T2D is characterized by high levels of glucose concentration in the blood. This increase in glucose levels can cause cardiovascular diseases and, if left untreated, will lead to organ failures. For T2D patients, low sensitivity to insulin, which is the hormone responsible for lowering glucose concentration in the blood, causes the beta cells in the pancreas to produce more of it to compensate. This will eventually weaken the cells and damage them which in turns will make the body fail to regulate glucose concentration [2]. Insulin based treatment is initiated at later stages of the T2D disease when changes in diets and physical activities accompanied with oral medications have failed. Clinically, it is difficult to calculate suitable insulin doses for each specific patient. Therefore, many patients experience uncontrolled hyperglycemia for a long period of time until they reach a safe level of glucose [28]. This happens due to the high variability of each patient and the different effect of their diet and lifestyle on insulin sensitivity. Furthermore, many patients suffer from hypoglycemic episodes if they overdose on insulin due to inaccurate calculations. Developing a model for the glucose-insulin dynamics helps with developing and testing algorithmic insulin dose guiders that handle patients variability better and ensure a safe reach for a desired blood glucose concentration. Moreover, having a model can helps in the development of hypoglycemic episodes alert systems for the patients. Additionally, models for glucose-insulin dynamics can help medical professionals with their patients' treatment plans.
In this paper, a survey of current models and simulators for the insulin-glucose dynamics in T2D patients is carried out. After that, a new model is proposed by combining and modifying some of the models from current literature. Finally, simulation results are provided to discuss the introduced model.

Literature Review
In general, there are two main categories of methods to model systems: first principles methods, or data driven methods derived by fitting data to general mathematical structures such as ARMAX models.
The glucose-insulin dynamical models for T2D patients based on first principles can vary with different degrees of complexity. Generally in the literature, there exist two main categories of such models: minimal models and maximal models [6]. Maximal models are very detailed models which model metabolic functions at a molecular level. On the other hand, minimal models are less detailed and rely mostly on compartments and mass balance equations. While maximal models provide a great level of accuracy, the amount of different data which is required to estimate parameters for the models is large and difficult to obtain from patients undergoing typical treatment plans. Moreover, the high accuracy of maximal models provides little relevance to the accuracy of the general glucose-insulin dynamics within the human body [6].
In contrast, minimal models consist of compartments to represent the distribution, diffusion, and production of glucose and insulin in the body with terms to represent the interaction between them. Furthermore, these models include pharmacokinetic equations to describe exogenous insulin injections and the intake of other medications. Several simulation oriented models for type 1 diabetic (T1D) patients were developed such as the ones in [19][16] [14]. The models for T1D patients can be extended to model T2D patients by including dynamics for endogenous insulin production and by obtaining new probability distributions for their parameters with the use of T2D patients data. As for already existing T2D models, Cobelli's model in [7] is a detailed model for glucose-insulin dynamics which is intended for simulation use. It includes nonlinear terms with hybrid differential equations. The model parameters are mainly estimated from people who do not suffer from diabetes. However, general parameters for T2D patients were also estimated using data from T2D patients. Although the model does not consider molecular level dynamics, its level of complexity still requires a variety of measured variables to estimate parameters for a given patient. In addition to glucose measurement in the plasma, the model used traced consumed glucose data, plasma insulin concentration measurements, and data regarding C-peptides (Amino Acids molecules which are side product of insulin secretion). The authors in [7] provided only mean parameters for T2D patients. They later developed an official simulator based on an extended version of their model to account for oral medications, Glucagon, and physical activities. The simulator is reported to have joint probability distribution for the parameters of T2D patients [27]. Nevertheless, these distributions were not published. Moreover, there are no published material for the equations and the structure of the extended model. Another detailed model is recently developed by [9]. It is an extension from the model in [23]. The model includes the affect of oral medications, Glucagon, and Glucagon-like peptide-1. The model is provided with mean parameters. Some of the parameters are identified from clinical data, while others are taken from different sources of literature. None of the mentioned T2D models so far considers or provide a mathematical structure for injected insulin, ingestion of multiple meals, and physical activity. A simpler T2D simulation oriented model with insulin injection is provided in [22]. It was developed as an extension to the model in [15] to take into account when both fast acting and long acting insulin are used for treatment. The model also uses modulator functions to account for the circadian rhythm in the glucose-insulin dynamics. As for parameter estimation, the model requires data of consumed glucose, glucose concentration in the plasma, injected insulin, and endogenous insulin concentration to be identifiable. In [22], parameter estimation was carried out using data from two different trials with different types of insulin with a total of 29 T2D patients. The parameter estimation took into account interindividual variety by estimating a mean and a variance for the parameters. Therefore, the model can be used to simulate a population of T2D patients. The model, however, is simple and does not include oral medications, physical activity, and glucose ingestion. Another model is the one from [1]. This model requires plasma glucose measurement with injected insulin concentration data. Therefore, this model can be validated, improved, and fitted with data collected from patients during typical treatment plans as done by the work in []. The model takes the T1D model in [16] and extended it with a term to describe endogenous insulin production. The model, nevertheless, does not include a mathematical structure for physical activity or oral medications. Probability distribution for the model parameters can be found in [1].
For data driven models, several methods were attempted for T1D in [13,10,26] and T2D in [20]. These models were generally developed to predict and detect hypoglycemic risks. However, data driven models are difficult for the purpose of developing a general model for T2D patients simulations. One major problem is that when data is collected from patients, the patients are already following some form of a feedback control mechanism between measured glucose levels and injected insulin. Therefore, fitting models on data from insulin inputs to glucose output will also include the dynamics of the feedback mechanism. This is undesired when the model is intended to test different control algorithms or if it is intended for the development of control strategies. This problem can be solved by perturbing the input insulin doses with time varying functions. However, it is difficult to impose perturbation on the insulin input for the patients without affecting their comfort and health.
In this work, it is intended to provide a model with mathematical structures for the effect of multiple glucose meals, insulin injections, multiple oral doses of metformin with different sizes, physical activity, and stress. The model is based on the one in [9] with modifications and inclusions as following (see figure 1): • Modifying the model to account for multiple meals (see section 3.1). • Including a model for insulin injections based on [18] (see section 3.2).
• Modifying the metformin model to account for multiple different doses in section 3.3. • Including the effect of physical activity based on [4] (see section 3.4). • Including the effect of stress based on [11] (see section 3.5).
In addition, two simulation cases are provided in section 4 to demonstrate the effect of lifestyle changes with the model. The model proposed in this paper is intended to be used as a starting point for developing an open source simulator for T2D patients when data is available. All the simulations in the paper are done with a Matlab code which can be obtained upon request from the authors. The model parameters which are used in the simulations are found in B.1.

Model Description
The model is mainly based on the one from [9] with the following four main subsystems: • Glucose subsystem.
See figure 1 for an overview of the model. The glucose and insulin subsystems are modelled as a set of compartments representing different main parts of the human body: brain, heart and lungs, guts, liver, kidney, and peripherals. The flow between these compartments follows the human blood cycle. As for the glucagon and the incretins, a single compartment is used for each one of them as it is assumed that glucagon and incretins have equal concentration in all the body parts. In addition, the model contains metabolic production and uptake rates for different compartments. These metabolic rates are generally defined as their basal values multiplied with scaling variables that depend on the concentrations of insulin, glucose, and/or glucagon (see (A.1)). The pancreas has a different nonlinear and hybrid model. In addition, a glucose ingestion model based on [7] is included as in [25] but modified to handle multiple meals along the day. Moreover, metformin and vildagliptin oral treatment models are included based on [24] and [17] respectively as in [9] but with a modification on the oral metformin model to handle different oral doses along the treatment. Additionally, a physical activity model based on [4] is added to the model. Furthermore, long acting and fast acting insulin injection models based on [18] are added. Finally, stress is included as a factor α s ∈ [0, 1] as in [11]. The main model includes parameters that were estimated by [23] for a healthy 70 kg male. The work in [25] considered a subset of these parameters to be estimated for the diabetic cases. Parameters for the different added models are taken from their corresponding literature. In the following subsections, the added and modified models and states will be discussed. The full model equations are provided in appendix A.

Glucose Absorption Model
In this section, a modification is introduced to account for multiple glucose meal sizes. The model used for glucose absorption in [9] considers only one glucose meal and was used for oral glucose tests where the patient were given an oral glucose dose and asked to fast while data is collected. The model is given as: Where q Ss (0) = D q [mg] is the oral glucose quantity, Ra is the rate of glucose appearance in the blood, f q is an absorption factor, k 12 [min −1 ] and k abs [min −1 ] are the rate constants for glucose transfer to stomach and glucose absorption in the intestines respectively, k empt [min −1 ] is a rate parameter for emptying the stomach of glucose to the intestines. This parameter can have values between k min and k max depending on the glucose dose size D q . In order to make the model handle different meals with different time instants, the parameter D q needs to be modified according to the meal sizes and time. The following are the proposed modifications:

Venous blood Arterial blood
G H , I H Heart and Lungs

Insulin Injection
Phyisical activity

Flow between compartments
Effect of a variable r Inj is the integer number of meals until time t, u qi [mg] is the amount of oral carbohydrates intake for meal i. The state D e is introduced to handle the accumulation of carbohydrates meals with a decay factor k q [min −1 ] in order to remove the effect of meals with time. With that, parameter D q is now a state updated by D e each time a new meal is consumed and made to converge to the last given meal amount u qi with the same rate factor k q such that it converges to the original model through time if no meal is consumed afterwards. Note that D q (0) = 0 when a zero carbohydrates meal (u q Nq (0) = 0) is assumed at time t = 0, which leads to (1d) being undefined (ϕ 1 q Ss = ∞ 0). To avoid this, the state D m (0) is set to 1 when u q Nq (0) = 0. Note that the value D m (0) can have any nonzero value in the case of zero carbohydrates meal at t = 0. The value D m (0) will not affect the rate of glucose appearance in the plasma since the states q Ss for ingested carbohydrates, and D e for the effect of accumulation of meals depend on u q Nq (0) and not D m . Parameters f q , k ϕ1 , and k ϕ2 are known and taken from [7]. The rest of the parameters, k 12q , k min , k max , k abs , and k q are to taken to be the mean parameters which were estimated in [25]. The introduced parameter k q has no estimate. Therefore, it is assumed to be equal to k min .
A simulation of a patient with the modified meals model compared against the unmodified one is shown in Figure 2.
for the insulin concentration in the interstitial fluid periphery compartment. It can be seen from the simulation results that the glucose appearance in plasma is distributed in a larger window of time with lower peaks for meals that are close to each other. This is due to the reduction of the stomach emptying rate k empt in response to increased accumulation of ingested carbohydrates captured by the state D q . Additionally, glucose appearance in plasma for the modified model in response to meals after hours of fasting closely resembles the glucose appearance in plasma for the unmodified model as can be seen for the dinner and breakfast meal. This is intended since the unmodified model was proposed for fasting conditions.

Insulin Injection Model
In this section, a model for long acting and fast acting insulin injections based on the one from [18] is introduced in [9]. Both fast and long acting insulin analogues treatments are considered for the model. When analogue insulin is injected, it dissociates from its hexameric form to dimers and monomers which then can penetrate the capillary membrane and get absorbed into the plasma. For fast acting insulin, only two compartments are considered: a compartment for insulin in its hexameric form, and a compartment for insulin in its dimeric and monomeric form. The following are the equations for fast acting insulin: Where r la , r f a ≤ 1 are the fractions of long acting and fast acting insulin that get to the periphery compartment, and V I P H [L] is the volume of the interstitial compartment. Figure 3 shows a simulation for a patient having the same basal values and following the same meal plan as the simulation discussed in section 3.1. The patient takes a long acting insulin dose of 50 [U] everyday an hour before the breakfast meal. Additionally, the patient takes a 30 [U] of fast acting insulin 15 minutes before dinner. Long acting insulin lower the glucose concentration over a large window of time. Moreover, the fast acting insulin helps at reducing the glucose peak after dinner.

Metformin
In this section, a modification for the metformin model in [9] is carried out to account for multiple doses of oral metfromin with different amounts. The metformin model used in [9], including the pharmacokinetic and its interaction with full glucose-insulin dynamical models, is based on [24]. The pharmacokinetic model of metformin in [24] is given as following: is the metformin amount in the gastrointestina lumen, parameters k go , k gg , k pg , k gl , k pl , k lp , k po [min −1 ] are transfer rate constants between the compartments, and M O is the flow rate of orally ingested metformin which is modelled as: Where A, B [µg min −1 ] and α M , β M [min −1 ] are constant parameters that were identified in [24]. These parameters were identified with data in which patients were taking only a 500 [mg] oral dose of metformin. Therefore, the model is modified in this work to take into account different amount of doses at different times by introducing the following:

Physical Activity Model
In this section, a physical activity model based on [4] is added to the model in [9]. The model in [4] was developed for a T1D model based on [3]. The model considers the change of the heart beat following a physical activity to be the stimulus of two states E 1 and E 2 which are dimensionless: Where t HR , τ e [min] are time constants, HR, HR b [bpm] are the current and rest heart rates respectively, and the parameters a e , n e are dimensionless parameters. The first state E 1 is used directly as a stimulus to increase the insulin-independent glucose uptake in response to a physical activity while the state E 2 is used for the longer lasting change of insulin action on glucose. The glucose and insulin model structure in [4] is simpler than the one considered in this work. Nevertheless, the inclusion of the physical activity for the model in this work is similar to how other models include physical activity, e.g., see [8]. With that, the effect of the state E 1 is included as an increase in the clearance rate of glucose in the periphery interstitial fluid compartment with a constant parameter β e [bpm −1 ] as 1 is the clearance rate for glucose in the periphery interstitial fluid compartment in (A.3h). The effect of E 1 can be removed to obtain the original model by setting β e = 0. As for the effect on insulin action, the state E 2 is introduced on the glucose metabolic rates which depend on insulin as following: • An increase in the periphery glucose uptake rate r P GU in the the interstitial fluid periphery compartment (A.3h) by a constant α e as (1 + α e E 2 ) r P GU . • An increase in the hepatic glucose uptake rate r HGU in the liver compartment (A.3e) with a constant α e as (1 + α e E 2 ) r HGU . • A decrease in the hepatic glucose production rate r P GH in the liver compartment (A.3e) with a constant α e as (1 + α e E 2 ) r HGP .
The effect of E 2 can be removed to obtain the original model by setting α e = 0. The parameters for the physical activity model are taken from [4] except for α e and β e which were tuned to have a similar effect to the ones demonstrated in [4,8]. Figure 5 shows a simulation for the patient described in section 3.1 when the patient exercise everyday before dinner raising the heart rate from a rest value of 80 [bpm] to a value of 140 [bpm] for 30 minutes. The immediate effect of physical activity is seen in the simulation results. In addition, the prolonged effect of physical activity on insulin action on glucose is seen in the figure.

Stress Effect
In this section, the effect of stress is included in the model [9]. In [11], the effect of stress was included as a multiplicative factor 1 + α s , with α s ∈ [0, 1], to the glucose and glucagon production rates. This is based  on the fact that stress causes a direct increase in the pancreatic glucagon production through an increase in catecholamines which in turns drives an increase of glucose production in the liver [29]. In addition, stress was also included as a multiplicative factor 1 − α s to the pancreatic insulin secretion rate based on [29]. Similarly in this work, the effect of stress is included in the model as following: • An increase in the plasma glucagon release rate r P ΓR in the glucagon compartment (A.9a) as (1 + α s ) r P ΓR . • An increase in the hepatic glucose production rate r HGP in the glucose liver compartment (A.3e) as (1 + α s ) r HGP . • A decrease in the pancreatic insulin release rate r P IR in the insulin liver compartment (A.10d) as (1 − α s ) r P IR Figure 6 shows the effect of stress in a simulation for the same patient discussed in section 3.1 when the patient is stressed on the second day with α s ramping up from 0 to 0.4 in 6 hours, staying at 0.4 for 12 hours, and then ramping down to 0 for the rest of the day. Stress manages to increase glucose concentration together with a decrease in insulin concentration.

Simulation Results
In this section, two simulation cases are carried out and discussed briefly to show how the model incorporates the effect of lifestyle changes on type 2 diabetic patients. The first case is for a patient having the same basal values and meal plan as the one discussed in section 3.1. Additionally, the patient experiences stress on the second day in the same manner as for the patient discussed in section 3.  plasma glucose concentration of 6 [mmol L −1 ] by the end of the third day, the patient in the second case has lower plasma glucose concentrations during the day when compared to the first patient.

Conclusion and Future Work
The proposed model is shown to be able to be used with multiple glucose meal sizes, different metformin doses, physical activity, insuling injections, and stress. The model, however, need to be confirmed with real patients data. Moreover, real patient data can be used to estimate joint probability distribution to model a population of T2D patients.

A Full Model Equations
The compartments include metabolic production rates r CXP and metabolic uptake rates r CXU for substance X in compartment C generally defined as following: Where M I ,M G , and M Γ are multiplicative quantities for the effect of insulin I, glucose G, and glucagon Γ respectively, and r b CXP,U is the basal metabolic rate of substance X in compartment C. The general form for the multiplicative quantities representing the effect of a substance Y in compartment C with concentration Y C is given as: Where Y b C is the basal concentration of substance Y in compartment C, and a, b, c, and d are model parameters.

A.1 Glucose Sub-Model
Applying mass balance equations over the compartments for glucose, the following equations are obtained: is the transcapillary diffusion time for compartment i, and r xP , r xU are metabolic glucose production and uptake rates respectively. The following are the meanings of each subscript in the model: , and peripherals M p [µg] respectively. These coefficients increase (or decrease) the glucose uptake (or production) as seen in (A.6). The equations for these coefficients are given as following: Where ν GW,max , ν L,max , ν P,max are parameters to represent the maximum effect of metformin in each one of its corresponding compartments, ϕ GW,50 , ϕ GI,50 , ϕ GI,50 [µg] are the masses of metformin within the different compartments to produce half of its maximum effect, and n GW , n L , and n P are shape factors.

A.2 Incretins Sub-Model
The incretins hormones are metabolic hormones released after eating a meal to stimulate a decrease in blood glucose levels. For T2D patients, Glucagon-Like-Peptide-1 (GLP-1) is the most active incretin [12]. GLP-1 is then modelled with the following two compartments as in [9]:

A.3 Glucagon Sub-Model
The glucagon subsystem consists of one compartment as it is assumed to have the same concentration over all the body: Where r P ΓR is the plasma glucagon release rate. The state Γ represent a normalized glucagon state with respect to its basal value. This is done since it is difficult in practice to obtain glucagon measurements for each subject in order to initialize the state. Therefore for this model, the basal glucagon state is 1.

A.4 Insulin Sub-Model
Applying mass balance equations over the insulin compartments will yield the following: Where r LIC , r KIC , and r P IC are the liver, kidney, and peripherals insulin clearance rates respectively and are defined as following: the pancreas insulin release is calculated by the following: Where S [U min −1 ] is the pancreas secreted insuln rate, and S b , r b P IR are the basal values. The model for S and S b is described in subsection A.5.

A.5 Pancreas Sub-Model
The model consists of two main compartments: a large insulin storage compartment m s [µg] and a small labile insulin compartment m l [µg]. The flow of insulin from the storage compartment to the labile insulin compartment is dependent on a dimensionless factor P with a proportionality constant γ [µg min −1 ]. The factor P depends on a dimensionless glucoseenhanced excitation factor represented by X and GLP-1 through a linear compartment with constant first order rate α [min −1 ]. Upon a glucose stimulus, the glucose-enhanced excitation factor X will increase instantaneously depending on the glucose increase in the plasma G H . In addition, a dimensionless inhibitor R for X will increase in response to X through a linear compartment with a first order constant rate β [min −1 ]. During that increase, the secreted insulin S will depend directly on both X and its inhibitor R together with GLP-1. Afterwards when R reaches X or X starts decreasing after R reaching it, the insulin secretion rate will only depend on X and GLP-1. The following are the equations of the model: Where m l0 is the labile insulin concentration at zero glucose concentration. This parameter in [23] is provided with a value of 6.33 [U] for a healthy 70 kg male.

A.6 Vildagliptin
The vildagliptin model is based on [17]. The absorption of orally ingested vildagliptin is modelled by two compartments as following: