Mathematical model of blood glucose dynamics by emulating the pathophysiology of glucose metabolism in type 2 diabetes mellitus.

Mathematical modelling has established itself as a theoretical tool to understand fundamental aspects of a variety of medical-biological phenomena. The predictive power of mathematical models on some chronic conditions has been helpful in its proper prevention, diagnosis, and treatment. Such is the case of the modelling of glycaemic dynamics in type 2 diabetes mellitus (T2DM), whose physiology-based mathematical models have captured the metabolic abnormalities of this disease. Through a physiology-based pharmacokinetic-pharmacodynamic approach, this work addresses a mathematical model whose structure starts from a model of blood glucose dynamics in healthy humans. This proposal is capable of emulating the pathophysiology of T2DM metabolism, including the effect of gastric emptying and insulin enhancing effect due to incretin hormones. The incorporation of these effects lies in the implemented methodology since the mathematical functions that represent metabolic rates, with a relevant contribution to hyperglycaemia, are adjusting individually to the clinical data of patients with T2DM. Numerically, the resulting model successfully simulates a scheduled graded intravenous glucose test and oral glucose tolerance tests at different doses. The comparison between simulations and clinical data shows an acceptable description of the blood glucose dynamics in T2DM. It opens the possibility of using this model to develop model-based controllers for the regulation of blood glucose in T2DM.

Introduction

For some decades, mathematical models have been used in biological sciences to understand diverse aspects of diabetes mellitus (DM)[1]. For example, DM progression[2,3], diagnostic test evaluations[4,5], long-term micro and macrovascular complications[6,7], and blood glucose dynamics[8-10], among others, have been modeled. Particularly, mathematical models to emulate blood glucose dynamics in DM have been classified, according to the complexity of their description, in two major groups[11]. The first group considers the whole-body models developed under a pharmacokinetic–pharmacodynamic (PKPD) approach, which is characterized by being structurally simple with a limited physiological interpretation. The second group considers the physiological based PKPD (PB-PKPD) models, which mathematically describe the physiological interactions between different subsystems of the human body. Due to its structural simplicity, most of the models in the literature are PKPD[1]. Although these models are widely used, they do not include most of the processes responsible for glucose homeostasis. Hence, its use to model complex processes in DM, such as the DM pathophysiology, is limited and, it induces a trend toward the development of PB-PKPD models[1]. These models have focused on emulating the metabolic processes involved in glucose homeostasis, and are usually organ-based. Moreover, the PB-PKPD models of blood glucose dynamics in type 1 DM (T1DM) have been useful to synthesize model-based controllers for blood glucose regulation in T1DM[12-16]. However, type 2 DM (T2DM) affects multiple subsystems of the body and, consequently, the mathematical representation of the metabolic abnormalities in T2DM is challenging[17].

One of the most widely used PB-PKPD models was performed by Sorensen[10]. This organ-based compartmental model emulates the blood glucose dynamics of a healthy human by considering the main glucose metabolic rates as mathematical functions. In this model, each mathematical function was individually fitted to a set of clinical data of healthy people where the metabolic response of the patients was measured for different stimuli. Then, the physiology of the main metabolic rates of a healthy human body was mathematically reproduced. Although Sorensen’s model is quite robust, it has some limitations. For instance, the model does not include blood glucose and insulin dynamics in the pancreas. Instead, a single function representing the pancreatic insulin release rate is connected to the bloodstream. The above does not represent the physiology of the human body. In addition, the model does not consider the effect of gastric emptying. Therefore, the blood glucose dynamics after oral glucose intake, and the potentiating-insulin effect of the incretin hormones cannot be reproduced.

An extension of the Sorensen’s model, which covers its main limitations, was proposed by Alverhag and Martin[9]. Thus, the model included two ordinary differential equations (ODE) to quantify, through mass balance, the time-variation of the blood glucose and insulin in the pancreas. Additionally, the gastric emptying process, and the enhancing effect on insulin due to the incretin hormones were included by considering two new subsystem attached to the model. Furthermore, Alverhag and Martin hypothesize that a model of the blood glucose dynamics in T2DM can be developed by identifying the parameters of the mathematical functions representing the metabolic rates related to the pathophysiology of this condition[9]. Based on the above, Vahidi et al. used a nonlinear optimization approach to identify some parameters of the Sorensen’s model from a single data set of an oral glucose tolerance test (OGTT) in T2DM patients[8]. Even though in this article, the system response acceptably reproduces the OGTT, the set of identified parameters that minimize the error between clinical data and the system may not be unique. Therefore, it cannot be assured that the metabolic functions containing the identified parameters emulate the pathophysiology of the T2DM individually.

Consequently, this article proposes a PB-PKPD model of the blood glucose dynamics in T2DM, where some mathematical functions representing metabolic rates of the body, are individually fit to emulate the pathophysiology of the T2DM. Moreover, the effect of the gastric emptying, and the enhancing effect of insulin due to the incretin hormones are included to reproduce the blood glucose dynamics after oral glucose intake. To achieve this, the mathematical model of the blood glucose dynamics in a healthy human body, proposed in Alvehag and Martin[9], will be described as a set of 28-dimensional ODE. From the ODE set, the mathematical functions representing the impaired metabolic rates in T2DM were individually fitted to clinical data of T2DM patients by using the least-squares method (LSM). The clinical data were taken from several clinical tests where direct measurements of the tissues or organ response to local changes in solutes concentration were made. The resulting model was numerically simulated to test its ability to reproduce the blood glucose dynamic in T2DM patients for different inputs, and initial conditions. Finally, the error between the simulation, and the clinical data of the T2DM patients is quantified by using a statistical function.

This manuscript is organized as follows: “Methodology" shows the methodology, while the results and discussions are set out in “Results and discussion”. Finally, the article ends with some concluding observations on this work in “Concluding remarks”. Also, Supplementary information have been included to show the nomenclature of all the variables and the numerical values of the parameters contained in the equations used throughout the manuscript.

Methodology

The mathematical model in Alverhag y Martin is a nonlinear dynamic system consisting of four clustered subsystems[9]. The subsystems are compartmental representations of the human body, where each compartment represents an organ or tissue where an important process of mass exchange is carried out. The compartments are interconnected through the blood flow. Then, by means of a mass balance in the compartments, each of the subsystems quantifies the concentration of one solute (i.e., glucose, insulin, glucagon, or incretins). A detailed explanation of the system and its nomenclature can be found the Supplementary information.

The system is a set of 28 ODEs composed of nonlinear continuous functions. Therefore, it follows that the solution of the system (x(t)) exists in a domain as long as the initial conditions are in . As a methodological approach in this work, the solution of the system is represented from a state-space theory as the vector:where is semidefined positive, which means that it belongs to the set . Using the state definition in Eq. (1), the system is defined as:where the vector field determines the time evolution of x(t) starting at initial condition () in the initial time (, and contains the parameters in the functions representing hemodynamical processes, while H contains the parameters in the functions representing the metabolic rates of the system. The parameter values of the system in Eq. (2) can be found in the Supplementary information.

Model simulation and initialization

The mathematical model in Eq. (2) successfully simulates the blood glucose dynamics of a healthy human body after intravenous glucose infusion and oral glucose intake[9]. For the above, an input to the system is considered containing: (i) a continuous intravenous glucose infusion rate (), which is introduced to the system as an insulin rate in mg(dLmin), and (ii) an oral glucose intake (), which is introduced to the system in mg and it is connected to the gastric emptying process (see the Supplementary information). The output of the system (y) is considered as and , whose meaning concerns to glucose and insulin vascular concentration in peripheral tissues, respectively. The time evolution of y is used to compare the model simulation with clinical data where the glucose and insulin concentrations are taken from a blood sample of the patient’s forearm during a test. For all the simulations, the model in Eq. (2) was numerically solved by using a variable step in the function ode45 (Dormand-Prince) of MATLAB[18]. The simulation time was defined as the time length of the clinical trial.

For model initialization, the basal condition and were computed from the solute concentrations in the fasting state of the patients. The condition is determined as the mean fasting glucose and insulin concentration from the blood samples collected over several days, this is and , respectively. The condition is determined as the fasting glucose and insulin concentrations from a blood sample at time zero of the clinical test; this is and , respectively. Mathematically, the fasting state has a physiological correspondence with the steady-state of the system () in Eq. (2), this is:then, since interstitial, arterial, and venous concentrations are the same at the steady-state, the peripheral vascular data for and are computed from the arterial or venous data. The remaining 26 components of and are obtained from the solution of the Eq. (3).

Metabolic rates of the model

The subsystems described in the Supplementary information are coupled by the functions representing the metabolic rates of the glucose, insulin, glucagon, and incretins. These metabolic rates are mathematically modeled as constant or linear functions of the mass accumulation in the compartments; or multiplicative functions of the metabolic basal rate. Specifically, the metabolic rates in the glucose, and glucagon subsystems are multiplicative functions with the following general form:where represents the basal value of the metabolic rate r, and each M is the isolated effect of the normalized concentration of glucose (), insulin (), and glucagon () of the normalized metabolic rate (). The above implies that when the glucose, insulin, and glucagon are basal, therefore . To represent the characteristic sigmoidal non-linearities of biological data correlations, excepting the isolated effects that are states of the system in Eq. (2) (i.e., and ), all the isolated effects are hyperbolic tangent functions of some normalized component of the state, this is:where for , and with are dimensionless parameters. A list containing the nominal values of the parameters can be found in the Supplementary information. Using these values, the system in Eq. (2) simulates the blood glucose dynamics after an intravenous glucose infusion or an oral glucose intake in a healthy human body[9]. For the mathematical modelling of the blood glucose dynamics of T2DM, the pathophysiology of T2DM must be emulated by modifying the value of the parameters of the functions representing the metabolic rates responsible of the characteristic hyperglycaemia. The above will be described in “Curve fitting”.

Curve fitting

For decades, different studies have identified the metabolic problems associated with the progression of T2DM in healthy humans[19,20]. It has been found that these problems are related to the metabolism of fats, and carbohydrates[19,20]. The metabolism of this latter is the object of study in this work.

Mainly, the pathophysiology of the T2DM is characterized by[19]: (i) insulin resistance, defined as an impaired effect of insulin on glucose uptake by peripheral tissues, (ii) excessive hepatic glucose production, due to accelerated gluconeogenesis, and (iii) -cell dysfunction, represented by an impaired pancreatic insulin release. Then, the mathematical functions of the system in Eq. (2) modelling the aforementioned metabolic rates are: the effect of insulin in peripheral glucose uptake (i.e., ), the effect of glucose, insulin, and glucagon on the hepatic glucose production (i.e., , and , respectively), and the pancreatic insulin release (i.e., ). Since a small variation in the parameters of the after-mentioned metabolic rates results in a variation of the solute concentrations in the model, in the following sections, the terminology of the sensitivity analysis from Khalil will be adopted[21]. Therefore, the above metabolic rates will be called sensitive metabolic rates.

In what it follows, the sensitive metabolic rates were selected to fit the clinical data of T2DM patients. Explicitly, the fitting of is supported by several clinical tests where a decrease of the first phase of pancreatic insulin release in patients with T2DM is exhibited[22-24]. The above is consistent with the early proposal to induce a partial impairment on insulin release from the labil compartment, in order to decrease the first phase of insulin release in T2DM patients[25]. Due to the above, the functions representing the first phase of insulin release (X and ), and the time-variation of the amount of labile insulin ready to be released, were studied by a sensitivity analysis as in Khalil[21] to select the parameters that show a major contribution to the sensitivity on solution . The selected parameters were identified from the clinical data of T2DM patients. The rest of the parameters remained unaltered.

Static and dynamic fitting approach

To solve the parameter fitting problem, two things are required: The set of clinical data used for the isolated effects fitting was obtained from selected clinical tests of T2DM patients. The conditions of each one of the selected articles are consistent with those originally considered for mathematical modeling in Ref.[10]. These conditions are compiled in Table 1. In the selected articles, the clinical data was taken from a set of individual with no other significant medical history than T2DM. Nevertheless, for curves fitting, we used the reported mean value of the tissue/organ response to local changes in the solute concentration of the subjects. Originally, to mathematically model the metabolic rate , Grodsky obtained data from a graded glucose step-response with the isolated perfused pancreas in rats[25]. Since it is impossible to obtain this data from humans, the selected parameters of this metabolic rate were identified using clinical data from an input–output approach of the system, in Eq. (2). The data were taken from an OGTT in DeFronzo et al.[26], where the plasma glucose and insulin response to oral intake were measured in nine T2DM subjects after the consumption of 1 g/kg-body weight of oral glucose.

Table 1

Conditions of the clinical test and its interpretation in the mathematical model.

Rate	Conditions
MPGUI	The glucose concentration was maintained in its basal state by a glucose clamp (i.e., MPGUG=1). The above allows to observe the isolated effect of insulin in the peripheral glucose uptake (i.e., MPGUI) by knowing the normalized value of the peripheral glucose uptake (i.e., rPGUN=r_PGU/rPGUB)
MHGPG	Somatostatin was administered to block the endogenous release of insulin, and glucagon. Exogenous insulin, and glucagon replacements were introduced to the patients to maintain them in their basal state (i.e., MHGPΓ=1,x_21=x21∗ ). The above allows to observe the isolated effect of glucose on the hepatic glucose production (i.e., MHGPG) by knowing the normalized value of hepatic glucose production (i.e., rHGPN=r_HGP/rHGPB)
MHGPI_∞	Glucose concentration was maintained at basal state by means of a glucose clamp. The above allows to observe the isolated effect of insulin in the hepatic glucose production after a stabilization time (i.e., MHGPI_∞)
MHGPΓ_0	Somatostatin was administered to block the endogenous release of insulin, and glucagon. Exogenous insulin, and glucose replacements were introduced to the patients to maintain them in their basal state by means of a glucose clamp (i.e., MHGPG=1,x_21=x21∗ ). The above allows to observe the isolated effect of glucagon on the hepatic glucose production (i.e., MHGPΓ_0) by knowing the normalized value of hepatic glucose production (i.e., rHGPN=r_HGP/rHGPB)

The clinical data from the studies that fulfill the criteria in the table were used to fit the isolated effects of the sensitive metabolic rates of the model.

A set of clinical data in T2DM patients.

A mathematical method to fit such data to the function representing the sensitive metabolic rates.

The mathematical method used to fit the functions to clinical data is the least squares (LSM). In general, the LSM lies that the following relation is fulfilled[27]:where z, and are vectors containing n observations, and is a vector of p unknown parameters of the sensitive metabolic rate. To estimate the n values of g are computed for all z. Then, is the estimation of the vector of parameters corresponding to that minimizes the residual sum of squares of an objective function over some feasible the vector of parameters . The isolated effects of the sensitive metabolic rates were fitted to clinical data by a static approach of the LSM. After that, a dynamical approach of the LMS was used to identify the parameters of the function. In what follows, both approaches will be described.

In the static approach, the unknown parameters from the Eq. (5) are grouped as . The vector is estimated with an iterative process using the following objective function:where is clinical data of the mean of the normalized metabolic rate in T2DM patients respect its basal value in Ref.[9], and is the clinical data of the mean of the normalized solute concentration taken from the forearm. The minimization of the objective function in Eq. (7) was numerically solved with the function lsqcurvefit of the optimization toolbox of MATLAB[18]. The iterative algorithm used to find was ‘trust-region reflective’ proposed in Li[28]. After fitting, (,) are graphically compared with the fitted isolated effects functions. Then, the values of the parameters in were replaced by the values in .

In the dynamical approach the selected parameters from were grouped as with . The vector was estimated with an iterative process using the following objective function:where , and are the clinical data obtained from the mean of glucose and insulin concentrations, respectively, taken at the time, the weights and are the mean of the basal glucose, and insulin concentrations, respectively; and , were obtained from the model simulation. The above clinical data was taken from DeFronzo et al.[26]. The LSM problem in Eq. (8) was numerically solved using the function fmincon of the optimization toolbox of MATLAB[18] with the iterative algorithm ‘interior-point’. After the identification of the parameters of , the values in (from the static, and dynamical approach) were replaced by in order to emulate the pathophysiology of T2DM. Hereinafter, the resulting model is called T2DM model.

Comparison of the T2DM model with clinical data

The T2DM model was numerically simulated for comparison with a clinical test in T2DM where the blood glucose dynamics is observed after different stimuli. Considering that the route of glucose entry into the body plays an essential role overall glucose homeostasis[26], the T2DM model was simulated for the following test: (i) a programmed graded intravenous glucose infusion test (PGIGI) to account for the rapid response of the intravenous infusions, and (ii) an OGTT considering a dose of 50 g of glucose (50 g-OGTT), and a dose of 75 g of glucose (75 g-OGTT) to account for blood glucose changes due to the gastric emptying process, and the effects of the incretin.

The clinical data used to compare the DMT2 model with a PGIGI test was obtained from Carperntier et al.[29]. In this test, the glucose was administered intravenously in a total of 7 subjects with DMT2 (i.e., ). Mathematically, this is that the glucose was supplied through while . The duration of the test was 270 min distributed as follows: a basal sampling period was considered were from 0 to 30 min, after this, the steps of intravenous glucose infusion were introduced as , and 8 mg (dL min) for a period of 40 min each one. The conditions for model simulation were mg dL, and mU .

The clinical data used to compare the DMT2 model with an OGTT was obtained from Firth et al.[30], and Mari et al.[31]. In these test, 50 and 75 g of oral glucose was consumed by a total of 13 and 46 subjects with DMT2, respectively (i.e., or ). Mathematically, this is that the glucose was supplied through while . For the OGTT, the duration of the simulation was 180 min. The conditions for model simulation were 50,000 mg, mg dL, and mU L, for the 50 g-OGTT. Further, the conditions for model simulation were 75,000 mg, mg dL, and mU L , for the 75 g-OGTT.

The difference between the clinical data, and the model simulation was quantified with the following statistical expression:where , and G is the glucose concentration taken from the T2DM patients at the time . All the clinical tests were different from those used for parameter fitting.

Declarations

The source of clinical data was obtained from publicly available sources, namely, recognized research journals and properly cited through the manuscript. No person was directly involved in this study as a source of clinical data.

Results and discussion

The clinical data that fulfill the conditions provided in Table 1 were taken from the references grouped in Table 2. The parameter set for each isolated effect of the sensitive metabolic rate can be seen in Table 3. Furthermore, in Fig. 1 it can be found a graphic representation of the curves that fit the isolated effects functions of the sensitive metabolic rates to the clinical data of the Table 1. As can be seen, the curves in Fig. 1 do not necessarily pass through the point . This is because the isolated effects of the metabolic rates were normalized with respect to the basal value of the metabolic rates in Alverhag and Martin[9], which correspond to a mathematical model of the blood glucose dynamics in a healthy human body. The above is justified by the fact that not all isolated effects of glucose, insulin, or glucagon on a metabolic rate are observed altered in T2DM patients. Since the metabolic rates are expressed as multiplier factors of the basal metabolic rate, the isolated effects that have not been observed altered in patients with T2DM will continue to be multiplier factors of the basal metabolic rate () of a healthy human body.

Table 2

References of the clinical studies.

Rate	References	n_p
MPGUI	DeFronzo et al.[26]	9
	Vaag et al.[32]	12
	Kelly and Mandarino[33]	15
	Capaldo et al.[34]	6
	Kalant et al.[35]	11
MHGPG	Hawkins et al.[36]	10
	Mevorach et al.[37]	9
	Nielsen et al.[47]	9
	Del Prato et al.[48]	9
MHGPI_∞	Staehr et al.[49]	10
	Groop et al.[38]	9
	Campbell et al.[39]	14
	Baron et al.[41]	10
	DeFronzo et al.[26]	9
	Revers et al.[40]	10
	DeFronzo et al.[50]	38
MHGPΓ_0	Matsuda et al.[42]	8
	Baron et al.[41]	10

The table shows the set of references containing the clinical data used to fit the isolated effects of the sensitive metabolic rates. Column indicates the number of patients analyzed in each reported clinical study according to the reference in the central column of the table. The proposed parametric adjustment results from taking the means of each set of patients.

Table 3

Vector of fitted parameters from the static approach.

M(z_k,θ^)	θ^
MPGUI(x15N,θ^)	[7.9869,7.2537,0.4852,-5.2518]
MHGPG(x4N,θ^)	[1.0720,-1.0064,0.8712,-1.4930]
MHGPI_∞(x12N,θ^)	[0.3240,-0.2020,0.7625,-3.6977]
MHGPΓ0(x18N,θ^)	[0,1.495,0.6773,-0.0469]

The table shows the parameter values that minimize the residual sum of squares of the objective function for the different isolated effects.

Figure 1

Isolated effects fitting to clinical data. In these plots the solid line represents the isolated effects functions (a) , (b) , (c) , and (d) fitted to the clinical data from T2DM patients. Each symbol represents the mean measured value of the tissue/organ response to a local change on the solute concentration, from subjects. For these metabolic rates, the fitting approach was static.

As can be seen in Fig. 1a the curve corresponding to goes close to the point (, . The above means that the insulin-stimulated peripheral glucose uptake in a T2DM patient does not differ much from the one in a healthy human when the fasting hyperglycaemia, and basal insulin concentration are maintained in the T2DM patient. This characteristic of the T2DM has been previously reported in several articles[32-35]. In contrast, as can be seen in Fig. 1b, for , the value of is higher than one. Considering basal hyperglycaemia, it means that the hepatic glucose production is higher in T2DM patients compared to that observed in healthy humans. The above has been previously reported by various articles where the effect of glucose on the hepatic glucose production rate was verified for healthy control subjects, and T2DM patients. Hawkins et al.[36] have associated this increase with accelerated gluconeogenesis since glycolysis is normal for healthy subjects and diabetic subjects. Besides, Mevorach et al.[37] report that this inefficient suppression is due to a deficient inhibition of glucose-6-phosphatase activity and/or lack of inhibition of glucose-6-phosphate formation.

The characteristic hepatic insulin resistance of the T2DM is evident in Fig. 1b,c. This can be observed in the behavior of the curves for high values of the solute concentration, where the hepatic glucose production can not be suppressed entirely despite significant increment of the normalized glucose and insulin concentration in the liver. The above is consistent with clinical evidence where the blood glucose has an impaired ability to inhibit the hepatic glucose production at basal insulin and glucagon concentrations in T2DM[19,36,37]; and the insulin concentration is ineffective to suppress the hepatic glucose production at basal glucose and glucagon concentrations in T2DM[26,38-40]. Finally, the role of glucagon in hepatic glucose production in T2DM patients can be seen in Fig. 1d. In this graphical representation, the behavior of the function is consistent with the clinical data of patients with T2DM[41,42]. Then, it follows that by fitting the isolated effect functions to the clinical data, it is possible to individually emulate the pathophysiology of the T2DM.

After isolated effects fitting, a parameter set containing the parameters of that shows a greater contribution to the sensitivity on the solution was selected. From the sensitive analysis, the selected set of parameters was . The values of that minimize the objective function in Eq. (8) can be seen in Table 4.

Table 4

Vector of identified parameters from the dynamical approach.

f(z_k,θ^)	θ^=[η_36,η_39,η_40,η_42,η_44,η_45]^T
x_6(t,θ^)	[3.2717, 2.8504, 0.9330, 0.0867, 7.6707, 0.0565]
x_14(t,θ^)

The table shows the parameter values that minimize the residual sum of squares of the objective function for the sensitive metabolic rate .

Once the nominal values of the parameters are replaced by those in of Tables 3, and 4, the system simulates the response to a 70 g-OGTT. The results of the simulation, and its comparison with clinical data from DeFronzo et al.[26] can be seen in Fig. 2. As can be seen there, is an acceptable approximation of the simulation curve to the clinical data. Moreover, during almost all the simulation time, the model output remains within the bars of the standard error. It should be noted that even when , and have different orders of magnitude, the emulation of both blood glucose and insulin was successfully achieved. This is mainly due to the addition of weight functions of weights in the objective function of the Eq. (8).

Figure 2

Graphical result of fitting to the clinical data. In these plots the solid line represents the variation of the (a) glucose or (b) insulin concentration in the peripheral compartment of the T2DM model. The symbols represent the mean±SEM value of the solute from the subjects. These data were taken from DeFronzo et al. where a 70 g-OGTT was performed[26]. For the simulation it was considered a consumption of 70 g of glucose at time equal to zero. The parameters are those whom minimized the objective function from the dynamic fitting approach.

As noted in “Methodology", since there is no clinical data of the individual response of measured against different stimuli, the dynamic approach used to fit is based on nonlinear optimization. As a result, the set of values obtained minimizes the objective function of the Eq. (8), nevertheless, it cannot be assured that the pathophysiology of the pancreatic insulin secretion in T2DM is individually emulated. However, due to the individual fitting of the isolated effects, the number of parameters to be identified by a dynamic approach is minimal. A proposal to avoid the above is to replace the pancreatic insulin subsystem with a model of the pancreas whose pathophysiology could be described by a set of clinical data of patients with DMT2.

After the metabolic rates fitting the resulting model (i.e., T2DM model) was simulated and compared with clinical data. In Fig. 3, it can be seen the T2DM model response for the PGIGI test, and the clinical data from Carpentier et al.[29]. As can be seen, the simulation of the T2DM model is not significantly different from the reported clinical data. Moreover, the absolute of the maximum difference between simulation, and the clinical data is 9.4 mg dL. This is consistent with the obtained statistical value mg dL for this test. It follows that the T2DM model can reproduce the step response of the blood glucose due to an intravenous glucose infusion input.

Figure 3

Simulation of a PGIGI test. In this plot the solid line represents the simulation of the blood glucose in the peripheral compartment of the T2DM model. The symbols represent the mean±SEM value of solute from the subjects. These data were taken from Carpentier et al. where a PGIGI test was performed[29].

Figure 4 shows the T2DM model response of the OGTT for different doses. Similarly to the observed clinical data, the model response rises to a maximum peak approximately at 80 min after the stimulus of 50 g and 75 g. The statistical value for the 50 and 75 g-OGTT is 16.84 mg dL, and 13 mg dL, respectively. As can be seen, after oral glucose intake, the response of the model for the OGTT test is relatively slow, showing a maximum peak at approximately 80 min after glucose stimulation. Compared with the results of the PGIGI, the increase in glycaemia in the OGTT is slower. This is because of the digestion process, after an oral glucose intake, induces a delay proportional to the glucose appearance rate in the gut. Furthermore, as can be seen in Figs. 2, 3 and 4 the basal blood glucose is slightly elevated compared to the concentration of a healthy subject.

Figure 4

Simulation of 25 g and 75 g-OGTT. In this plot the solid and dashed lines represent the simulation of the blood glucose in the peripheral compartment of the T2DM model during a 25 g and 75 g-OGTT, respectively. The symbols represent the mean value of solute from the n subjects. Particularly, the triangles represent the mean ±SEM. The data for the 25 g and 75 g-OGTT were taken from Mari et al., and Firth et al., respectively.

According to the World Health Organization guidance for diagnostic tests of DM, a fasting glucose concentration 126 mg dL is characteristic of DM[43]. Moreover, the patients with this impaired blood glucose should undergo by a formal 75 g-OGTT for DM diagnosis[43]. A representation of this test can be seen in Fig. 4 where after two-hour postload glucose the dotted curve had shown a glucose concentration 140 mg dL. This is a characteristic behavior of DMT2 that contrasts with that of a healthy subject, where the normal homeostatic glucose process results in a concentration of less than 140 mg dL after 2 h of the glucose intake.

Based on the values of the function, it can be concluded that the model emulates with acceptable precision what is reported in the clinical data for PGIGI, and OGTT. However, this model considers only the carbohydrate metabolism but not fat, and protein metabolism. Therefore, the effect of free fatty acids, and the physiology related to amino acids level on blood glucose dynamics are not included. Besides, the model does not consider the counter-regulatory effect of growth hormones, adrenaline, or cortisol. Nevertheless, the above can be considered later in the model by adding other subsystems for the free fatty acids dynamics, and other metabolic functions to consider the effect of the missing hormones.

Concluding remarks

The main contribution of this article was derivating a model for T2DM, including physiological features to emulate blood glucose dynamics. The modelling departs from a PB-PKPD modelling approach, and the individual fitting of the sensitive metabolic rates allows us to capture the pathophysiology of the metabolic rates in T2DM. This methodological procedure enables us to successfully emulate the blood glucose dynamics of T2DM after a continuous intravenous glucose infusion an oral glucose intake. As convincing numerical evidence of the above, Figs. 3 and 4 show to what extent the T2DM model predicts the clinical data.

The individual fitting of the sensitive metabolic rates to clinical data ensures that the pathophysiology of T2DM is preserved, such that diverse scenarios might be predicted. For instance, this model can be used to determine appropriate oral therapy for blood glucose regulation by connecting a PKPD model of a hypoglycaemic drug (e.g., sulfonylureas, biguanides, thiazolidinediones, among others). Such is the case of metformin therapy, where the target metabolic rates were modified by adding a multiplicative factor in , , and [44]. Similarly, the mathematical model we present lights up some complementarity to other research approaches. For example, the recent finding in improved glucose metabolism due to continued treatment with deuterium-depleted water (DDW) content in patients with T2DM[45] could be emulated with an organ-based model. Again, it is feasible to consider a multiplicative factor to the metabolic rate to reproduce the alteration on peripheral glucose disposal, as indicated in some clinical researches[46]. Furthermore, this model can be used to develop a feedback model-based controllers for blood glucose regulation in T2DM patients. This idea triggers the possibility to achieve the normoglycaemia by means of single or combined therapy of oral hypoglycaemic agents with an exogenous insulin input connected to . Finally, as a consequence of its mathematical structure, it is possible to consider structured or unstructured uncertainties in the described physiological-based model. Therefore, we can employ robust control techniques such as theory.

Supplementary Information.