Bayesian analysis of Glucose dynamics during the Oral Glucose Tolerance Test (OGTT)

1.   Introduction

This paper builds on the model of blood glucose dynamics during the Oral Glucose Tolerance Test (OGTT) of Kuschinski ^[1]. Our model considers the action and timing of insulin and glucagon in glucose homeostasis after an oral stimulus. We use the Bayesian paradigm to infer kinetic rates, namely insulin and glucagon secretion, gastrointestinal emptying, and basal glucose concentration in blood. We identify two insulin scores related to glucose concentration in both blood and the gastrointestinal tract, suggesting a classification for individuals with impaired insulin sensitivity.

Most tissues and organs in our body use glucose as an essential source of energy. The process of maintaining blood glucose at a steady state is called glucose homeostasis ^[2]. Food ingestion, fasting, physical activity, or exercise constantly shift the blood glucose concentration away from equilibrium throughout the day. A low blood concentration of glucose can causes seizures, loss of consciousness, or even death. On the other hand, the long-lasting elevation of blood glucose concentrations can result in blindness, renal failure, vascular disease, and neuropathy. Therefore, blood glucose concentration in healthy individuals stays within narrow limits. Two hormones secreted by the pancreas: insulin and glucagon, are the primary regulators of blood glucose. The World Health Organization (WHO) refers to diabetes type 2 as a chronic and metabolic disease characterized by elevated blood glucose levels. Type 2 diabetes is more common in adults and occurs when the body becomes resistant to insulin or does not produce enough insulin. One standard test to screen for diabetes is the Oral Glucose Tolerance Test (OGTT). In this procedure, the patient must drink a sugary solution after eight-hour of fasting. Several blood samples are taken over the process, the first before consuming the sugary solution, the rest are taken afterward to monitor the body's reaction to the glucose intake. A prediabetic state is characterized by an impaired fasting glucose level (IFG), an impaired blood glucose level at two hours after consuming the drink, referred to as impaired glucose tolerance (IGT), or both. Studies suggest that IFG may be associated with impaired insulin secretion while IGT is related to insulin resistance ^[3]. OGTT based diabetes diagnostics suffers from uncertainties along the test ^[4] on the amount of glucose absorbed and its absorption rate. A popular approach to describing glucose dynamics is compartmental models ^[5,6]. The inference of model parameters from glucose observations is called an inverse problem. The critical task of solving the inverse problem is to reliably predict quantities of interest in terms of the model´s variables and parameters. Of note, inverse problems are an imperfect path to knowledge due to errors in observation, modeling, and numerical simulation ^[7]. The Bayesian paradigm allows solving the inverse problem producing predictions with quantified uncertainties. Under this approach, the solution of the inverse problem is a probability distribution for the parameters of interest ^[8], which describes multiple consistent scenarios to fit the data. In this work, we use the Bayesian approach to model and analyze OGTT data. Here the regressor comes from an ODE system that describes the glucose dynamics. We model insulin and glucagon secretion originated from blood glucose concentration. Also, we consider delays and hormonal effects from the gastrointestinal tract.

1.1.   Related work

Research on glucose dynamics through compartmental models is quite common; see ^[5,6] for reviews. Usually, the description level is related to the crucial features of the diagnostic tests involved. For the Oral Glucose Tolerance Test, one of the main features to be modeled is the gastrointestinal tract ^[2,4,9]. Most models include the insulin effect on the glucose homeostasis process, include glucagon action is less common. Several works describe OGTT dynamics with a second-order differential equation ^{[10,11,12,13]}. Authors in ^[12] have shown that a linear model may describe some postprandial glucose excursions.The process in the gastrointestinal tract has become relevant in the last years. Incretins are hormones secreted by the gut, which stimulate insulin secretion before blood glucose level increases. This effect is reduced or even vanished for type 2 diabetes patients, see ^[14,15]. Also, gut hormone secretion is low in the fasting state ^[2]. Reactions of the body to the rise of the glucose level come with a delay, modeled by explicit delays in time ^[16,17] or by extra compartments in the dynamical system ^[1,4,18]. In epidemiology, latency and infection periods are usually described by an exponential distribution. Authors in ^[19] argued the inadequacy of this distribution since this proposal overestimates shorter and longer durations of the phenomenon in the period. A gamma distribution may describe more realistic distributions. We can obtain this effect by subdividing the compartment into $n$ stages. Specifically, the distribution obtained is an Erlang distribution, a gamma distribution with a shape integer parameter ^[20]. Inference problems look for determining the value of model parameters based on observed quantities of interest, like glucose, insulin, GLP-1, and/or GIP ^[4,18]. For a complex model, often, several sources of data are required to infer the parameters. Using just glucose data limits the complexity and may lead to an identifiability problem of the parameters. Continuous Glucose Monitoring (CGM) provides information about the body reaction during a more extended period ^[12,21,22]. Authors in ^[22] proposed a minimal model with CGM data and prove the structural identifiability using software packages. Authors in ^[11] proposed lifestyle adjustments for T2D subjects using glucose data for more than 4 hours each 30 min. The results show that the recovery time of post-prandial blood glucose level can be adjusted to 4 hours. Bayesian approaches for estimate parameters of compartmental models were proposed in ^[1,13]. Authors in ^[13] proposed a minimal model for glucose-insulin dynamics. They proposed interpreting the glucose-insulin system as a damped harmonic oscillator and suggested a classification based on two parameters: the peak of the curve and an average of rates. Authors in ^[1] proposed a Bayesian experimental design for time's sample collection. Since the inference results depend strongly on the quality of the data, an experimental design looks to improve the ability of the test to provide a more accurate result. Authors define a utility function for different time collections that allow determining an optimal among these candidates. For a Bayesian experimental design, this utility depends on posterior quantities.

1.2.   Contributions

The main contributions of this work are a model that includes: (ⅰ) A description of insulin secretion dependent on direct blood glucose level as well as gastrointestinal glucose. The incretin effect has a significant role in insulin secretion when the body reacts to an oral stimulus. Separately, these two sources allow us to propose a classification between patients. These parameters may suggest possible alterations in healthy patients and determine misclassified subjects. We obtain those by a Bayesian inference that takes less than five minutes and five glucose measurements, making it a fast, reliable, and accessible tool. (ⅱ) Delays of the body in the gastrointestinal tract and the endocrine system are modeled by gamma distributions. We introduce Erlang distributed periods in the hormonal dynamics for incretin, insulin, and glucagon.

1.3.   Limitations

We remark some limitations for this approach: (ⅰ) Since the glucagon dynamics turn on when the glucose level is lower than the basal level ($G_b$), for patients with all measurements higher than $G_b$, we do not have information about the parameter corresponding to this reaction. The consequence is that the posterior marginal for $\theta_2$ matches the prior marginal. Monitoring the glucose level for more than two hours may circumvent this situation. (ⅱ) The classification proposed is based on two parameters. Each is the product of two quantities with biological meaning. Nevertheless, we can not recover each action separately. Adding insulin data in the inference may allow recovery of each action parameter.

The organization of the rest of this paper is as follows. In Section 2.1, we present the ODE system that describes the biological situation for the glucose dynamics during OGTT. This section includes the description of the Bayesian inference problem in Subsection 2.2. Our numerical results are present in Section 3. Finally, we discuss our results and findings in Section 4.

2.   Materials and methods

2.1.   Modeling of Glucose dynamics during OGTT

In this section, we describe the physiological situation in detail to justify our approach. The biological modeling on which we rely is given by the following ODE system

where $G$ denotes the glucose blood level (in mg/dl), $I$ the insulin level above the basal insulin (in $\mu$U/ml), $L$ the glucagon level above the basal glucagon (in mg/dl), and $V$ the glucose from the sugary drink in the gastrointestinal tract (in mg/dl). We propose a dynamical system centered in the basal glucose level denoted by $G_b$, which typically takes values in the range $[75-100]$ mg/dl, ^[6]. Above $G_b$, the dynamics are regulated mainly by insulin action, while below $G_b$, mainly due to glucagon. This system describes linear interactions for both cases. A linear description for the load path through the gastrointestinal system is usual, while nonlinear functions describe glucose appearance in blood ^[18]. In this case, we decided to maintain the linearity assumption as in ^[10], which simplifies the biological phenomenon but reaches the objective of predicting the glucose level, a major quantity of interest for the inference ^[7]. Note that the complexity of models in ^[4,18,23] comes with a diversity of observed data to validate it. In our case, we proposed the inference using just glucose data at five times, which is a limitation to increase the complexity. Also, the model's simplicity is not a limitation to obtain a good fit and reliable inference ^[7]. Equation (2.1) models the glucose dynamics, which decreases proportionally to the insulin level and increases according to two sources, first by the glucose of the sugary drink from the intestinal tract, second by the glucagon action. $\lambda_1$ and $\lambda_2$ are efficacy rates for the glucagon and the insulin, $\lambda_3$ is the rate of glucose absorption from the gastrointestinal tract. The dynamics of insulin and glucagon are very similar and depend on the basal glucose level of the body. The pancreas secretes insulin when the blood glucose level is higher than the basal glucose level on the body. Conversely, the pancreas secretes glucagon when the blood glucose level is lower than the basal glucose level on the body. We model this secretion action by a switch on the quantity $G-G_b$. Mathematically speaking, we use the positive part of $G-G_b$, denoted by $(G-G_b)^+$, and defined by

to model the secretion of insulin, conversely $(G_b-G)^+$ to the secretion of glucagon. $\lambda_4$ and $\lambda_6$ are secretion rates for the insulin and the glucagon due to blood glucose level, and $\lambda_5$ and $\lambda_7$ are its disintegration rates. $\lambda_8$ is the insulin secretion rate due to the glucose level in the gastrointestinal tract. After a meal, insulin secretion may occur in two phases: an initial rapid release and a long-term release if glucose concentrations remain high ^[2]. Hormones from the gut play a main role in insulin secretion on oral glucose consumption. The second term on equation (2.2) mimics the secretion due to a high concentration of glucose in the gastrointestinal tract as for example, the incretin effect. The Incretin effect is reduced or even vanished in diabetic patients ^[14,15]. Therefore, insulin secretion may be due to direct glucose blood level and/or gastrointestinal glucose level. The model defined by equations (2.1)- (2.4) has eight parameters. Parameters $\lambda_5$ and $\lambda_7$ will be taken from the literature ^[1,25]. For the other six parameters, we will face an identifiability problem; that is, the structure of the model may lead to issues related to the uniquely estimation of parameters during the inference. To face this problem, we introduce scaled versions of insulin and glucagon. Note that by multiplying equation (2.2) by a constant, we obtain a scaled insulin

If we name $I_s = \lambda_I I$, we have $\dot{I_s} = \dot{\left(\lambda_I I\right)} = \lambda_I \dot{I}$. By substituying in equation (2.6), we obtain

We will apply this idea to obtain a scaled insulin and glucagon by multiplying equation (2.2) by $\lambda_2$ and equation (2.3) by $\lambda_1$. As we will explain later, we use only glucose observations for the inference, for which the loss of meaning for these scaled amounts do not represent a limitation. Finally we obtain the system

where $L_1 = \lambda_1 L, I_1 = \lambda_2 I, \theta_1 = \lambda_2 \lambda_4, \theta_3 = \lambda_2 \lambda_8$ and $\theta_2 = \lambda_1 \lambda_6$. Note that these substitutions allow us to decrease the number of parameters in the system, which will be helpful in the inference to fight the identifiability problem. Table 1 summarizes information about the model's parameters given by equations (2.8) - (2.11). For this model and the specific case when $G-G_b \geq 0$, we can show by the Similarity Transform Method ^[26] that $\theta_0$ is identifiable and $\theta_3$ is not, see supplementary section B. Note that, for the cases of OGTT curves with all data greater than $G_b$, the term $(G_b-G)^+ = 0$ during most all the test, which causes a non-identifiability of $\theta_2$ and $G_b$. Let us recall that the unidentifiability of parameters produces unrealistic values or high uncertainty around the estimated values of the parameters. A Bayesian approach deals with these problems by proposing a prior distribution on the parameters and approximating integrals from the MCMC, which produces less sensitivity on the results ^[27].

Also, we would like to address the possible delays of the body reaction. Introducing explicit time delays into the ODE system is very common ^[16,17]. The same effect may be obtained by introducing extra compartments. In epidemiology, these additional compartments model aspects as latency and infection periods with a more realistic distribution. This technique is known as an Erlang model ^{[19,20,28,29]}. Typically, people use the exponential distribution. Nevertheless, the exponential distribution overestimates the number of individuals whose duration of infection is shorter or longer than the mean ^[19]. This approach subdivides stage $I$ (and/or $E$) into identical substages in a classical SEIR compartment model. This modeling matches the renewal approach ^[20], which considers the cohorts of infectious(exposed) individuals while the ODE approach considers each infectious(exposed) individual. In this work, we consider the same idea to model hormonal behavior during the OGTT test. That is, instead of considering an exponential decay of quantities $V_1, L_1$ and $I_1$, we introduce extra compartments to obtain a different and delayed behavior of these into the bloodstream. The main difference is illustrated in Figure 1 (a) for the digestive compartment. An exponential decay ($m = 1$) gives higher weights to initial times. Note that the parameters $\lambda_5$ and $\lambda_7$ are the mean-life clearance of hormones insulin and glucagon. Under this approach, these parameters retain their meaning ^[20]. We propose a modification of the model in ^[30], deduced before in equations (2.8)-(2.11). This modification is given by

where $\theta_0 = \lambda_3$. After ingestion, glucose is absorbed in the upper gastrointestinal tract, transported to the liver, and finally reaches the peripheral circulation ^[23]. Introducing a second compartment for $V$ allows us to model this behavior of glucose due to the sugary drink in the digestive system as in ^[1,31]. Compartments $V_1$ and $V_2$ model the stomach and the small intestine. The second compartments for $I$ and $L$ allow us to model the delays due to the pancreas reaction ^[16]. These delays are a consequence of recurrent inhibitory dynamics ^[32,33]. In recurrent inhibition, we can see how the activation of a quantity produces excitation in a second quantity that inhibits the first's activity. This phenomenon is present in the dynamics of glucose-insulin-glucagon. The specific case for glucose-insulin dynamics is known as a negative feedback loop ^[34]. Following this approach, we can introduce extra compartments without new kinetic constants. This model is locally stable around the unique equilibrium point $(G_b, 0, 0, 0, 0, 0, 0)$ for subjects with small $\theta_1$. We show details about the stability of the model in Appendix C. In the following subsection, we propose a regressor from the dynamical system in equations (2.12)–(2.18) to fit the glucose data provided by the OGTT. The inference process follows a Bayesian approach.

2.2.   Bayesian Formulation for the Inference Problem

In this section, we consider the inverse problem of determining the posterior distribution for the parameter $\Theta = (\theta_0, \theta_1, \theta_2, G_b, \theta_3)$. Under the Bayesian approach, the solution of the inverse problem is a probability distribution conditioned on the information of the data. We have glucose measurements at five times $t_i = 0.0, 0.5, 1.0, 1.5$ and $2.0$ hours. We assume that the data $\textbf{y}$ follow the noise model

where $\mathcal{G}(\Theta)_i = G(t_i, \Theta)$ and $G$ is the solution for the glucose on the ODE system given in equation (2.12) and $\eta_i \sim \mathcal{N}(0, \sigma^2)$. The parameters $\lambda_5, \lambda_7$ are determined from previous works ^[1,31]. We propose $\sigma = 5$ mg/dl as in ^[1,13].

Our prior knowledge about the parameters are that $\theta_i > 0$ and the basal level of glucose is in the range of $[75,100]$ mg/dl ^[6]. We assumed independence on the prior parameters, that is

We propose gamma distribution for each parameter. Le us recall that if $Z \sim \Gamma (\alpha, \beta)$ then $\mathbb{E}[Z] = \alpha / \beta$ and $\mathbb{V}ar[Z] = \alpha /{\beta^2}$. Our priors proposal are given by:

The prior distribution for the parameter $\theta_0$ corresponding to the gastrointestinal dynamic is truncated. From simulations, we propose values greater than 0.5 to avoid almost constant trajectories for the glucose. For $G_b$, we consider a prior with mean $90$ mg/dl and variance $20$ mg/dl. The priors for $\theta_1, \theta_2$ and $\theta_3$ are less informative with same mean and variance in 10. Having the likelihood $\pi(\textbf{y} | \Theta)$ and the prior $\pi_0(\Theta)$, and since we are working in finite dimensions, Bayes theorem ensures the posterior distribution $\pi(\Theta|\textbf{y})$ existence and

Since we are using a numerical regressor in equation (2.19), the posterior distribution has no closed-form, and we need to generate samples of it. Note from equation (2.22) that the posterior distribution is known up to a constant of proportionality. Markov Chain Monte Carlo (MCMC) allows obtaining samples from a distribution under this condition. Let us recall that the chain is a sequence of random variables obtained from a transition kernel and distributed according to a stationary probability distribution. To illustrate the results, we consider estimators for this distribution. The more popular estimators are the conditional mean (CM), which is the mean of the posterior distribution; the conditional median, which is the median of the posterior distribution; or the maximum a posteriori (MAP), which is the value of the parameters where the posterior reaches its maximum value. These estimators are computed based on the samples generated from the MCMC. To generate samples, we used the t-walk ^[36]. We perform 10000 iterations for each patient with a burnin of 1000. To assess the efficiency of the sampler, we compute the integrated autocorrelation time (IAT). IAT is the number of steps required for the chain to produce independent samples. Authors in ^[37] argue that $IAT/n$, where $n$ is the number of parameters, is a measure of efficiency among MCMC samplers. Let us recall that the Monte Carlo error introduced in the approximations of integrals is proportional to $\sqrt{IAT/N}$, where $N$ is the number of samples from the posterior distribution ^[38]. That means we are looking for small values of the IAT to have efficiency. We show our results in the following subsection.

3.   Results

In this subsection, we show the results of the inference. We perform a MCMC using t-walk. Data were collected from 2012 to 2019 in the General Hospital of Mexico Dr. Eduardo Liceaga under the research Identificación de factores que predisponen a la diabetes mellitus tipo 2 en sujetos normoglicémicos con historia familiar de diabetes y su relación con obesidad. These data were previously used in ^[13]. All participants signed a consent. A first measurement is taken at fasting, after eight hours of fasting. After that, the patient has five minutes to consume a drink with 75 gr. of dextrose. Four different measurements of glucose are taken at 30, 60, 90, and 120 minutes. We have data from 80 female patients classified as follows:

1. 51 healthy patients (H).

2. 4 patients with Impaired Fasting Glucose (IFG): Fasting blood glucose level $\ge 100$.

3. 15 patients with Impaired Glucose Tolerance (IGT): Blood glucose level $\ge 140$ at $t = 2$.

4. 7 patients with IFG and IGT (IFG-IGT: both alterations).

5. 3 patients with Diabetes Mellitus 2 (T2D): Fasting blood glucose level $\ge 126$ and blood glucose level $\ge 200$ at $t = 2$.

Since only a fraction of the ingested glucose appears in the blood ^[39], we propose a value of $V_0 = 400$ based on the maximum value of the data set (375 mg/dl). The initial condition for the glucose level $G(0)$ is equal to the fasting observation $y_0$. Figure 1 (a) illustrates Erlang's distribution effects for $1, 2$, and $3$ compartments in the digestive tract. The case $m = 1$ corresponds to an exponential distribution. This case may describe the process of the glucose level in the stomach. For $m = 2$ and $3$, the quantity of glucose at $t = 0$ is 0 and may describe different parts of the intestine. We proposed the case $m = 2$ based on simulation results. Figure 1 (b) shows the RMSE at the MAP estimate for each category of patients in boxplots. The results are consistent with the standard deviation proposed on the likelihood ($\sigma = 5$). In Figure 1 (c)-(d), we show trace plots for the parameters $\theta_1$ and $\theta_3$ for several condition patients. These results illustrate the chain values without the burnin. In Figure 2, we show the fit to the data for several patients. We include (in grey) trajectories generated from posterior samples and the curve generated by the posterior median. The CM and the MAP estimates generate very similar trajectories.

From the inference information, we would like to explore the values of the parameters that describe different situations on patients. In Figure 3 (a), we show chain values for five different scenarios. In Figure 3 (b), we show the corresponding glucose measurements. Note that parameters $\theta_1$ and $\theta_3$ allow us to recognize low secretion insulin levels. The gastrointestinal emptying parameter $\theta_0$ allows us to identify the initial slope during the test. That is, $\theta_0$ determines the increase from level at fasting and level at $t = 30$ minutes. Two patients with low gastrointestinal emptying labeled as impaired2 and diabetic2 present higher values close to two hours instead of at the beginning. Note that this parameter does not allow us to recognize anomalies since, for two diabetic patients, the parameter may take very different values.

In Figure 5 (a)–(c), we show the posterior variance for the parameters $\theta_0, \theta_1$ and $\theta_3$ inferred from the MCMC. Note that for $\theta_1$, there are some healthy outliers with high variance, consequence of the unindentifiability problem, see details in Appendix A. In Figure 4 (a), we show the values of the MAP estimate for parameters $\theta_1$ and $\theta_3$ inferred. We plotted $1/\theta_3$ against $1/\theta_1$. A high value for $1/\theta_1$ may be produced by a high glucose level at time 30 minutes whilst a high value for $1/\theta_3$ may be produced by a high glucose level at time 1 hour and 30 minutes. Patients diagnosed with IFG are considered healthy in the classification based on insulin secretion indicators. Its condition is not a consequence a the oral stimulus. We performed a linear Support Vector Machine (SVM) algorithm using the package scikit-learn in python ^[40]. An SVM is a supervised machine learning algorithm for classification, which has become popular in biological applications ^[41]. In this case, we are looking for a separating line between two different classes, healthy patients and patients with an impaired condition. This hyperplane is optimal in the sense that it results from an optimization problem. Since naturally there is a transition stage from a healthy patient and a diabetic patient, we performed this algorithm for all the quantiles from 10 to 90 to include the uncertainty of our simulations in the classification. In Figure 4 (b), we show the glucose data, the fit, and the uncertainty corresponding to the patient with IGT inside the healthy zone in (a) (cyan star at left down corner). For this case, the parameter $G_b$ is identifiable, see Figure 5 (d). We show results of the inference for this patient in Table 2. Our results suggested that this subject may be misclassified as IGT instead of healthy.

4.   Discussion and conclusions

In this work, we propose a Bayesian approach to analyze OGTT data. The modeling includes two insulin indicators, one related to blood glucose level, the other with the glucose level in the gastrointestinal tract. These indicators describe insulin dynamics due to oral stimulus. Figure 4 (c) illustrates a possible classification for patients. Now, we discuss its scope and limitations. Patients diagnosed with IFG have an impaired value at fasting state and unimpaired values at two hours. Their insulin indicators classify these patients as healthy because their situation is not a consequence of the oral stimulus. A patient with IGT may have healthy values for its insulin indicators, as shown in Figure 4 (a) (cyan star in the left down corner). From its glucose values (Figure 4 (b)), we can deduce that its situation is a consequence of impaired glucagon secretion. Also, note that the glucose measurement at $t = 2$ is almost on the limit of a healthy condition. The insulin indicators allow recognizing this possible misclassification. Note that the detection of these cases is conditioned to have more data than fasting and 2h level sample. Finally, the classification places some healthy patients in the transition zone. This situation may result from a high glucose level at 30 minutes, for $\theta_1$, or at 1 hour and 30 minutes for $\theta_3$, or a measurement close to 140 at 2 hours. This situation may be the beginning of an anomaly. Control and follow-up of these patients for several months or years may confirm it. The scores presented here have physiological meaning, discriminate between healthy and diabetic patients, and are accessible since the inference takes less than five minutes and the data needed to perform it are just from glucose measurements. From our results, it seems to us that the classification in healthy, IFG, IGT, IFG-IGT, and T2D could be modified to recognize more specific details about the cause of an impaired condition. IGT subjects may suffer from diminished glucagon suppression and an impaired glucose level at 2 h may be caused by this irregularity ^[42]. Recognize the real causes of an impaired condition leads to different therapies and savings of health care resources. Also, early detection of anomalies and/or high risk of developing T2D allows for implementing healthy habits therapeutic strategies and postponing medication. The model described here is very simple but incorporates all essential elements of the physiological situation. As we mentioned before, the lack of complexity corresponds to the limited diversity of data. In this case, just with five glucose measurements, a very basic observation scenario, we were allowed to propose a classification. This classification enables recognizing possible misclassified patients and determining a transition zone, which may represent the risk of developing anomalies for healthy patients. Future work contemplates how to adapt the model to incorporate insulin and glucagon data. This information may provide possible explanations for an impaired condition.

Acknowledgments

We would like to thank J. Andrés Christen for his help and comments and Antonio Capella for valuable comments. Also, we would like to thank Adriana Monroy for the data and comments, and three anonymous reviewers for their comments to improve this manuscript. HF-A was partially found by RDCOMM grant from J. Andrés Christen. The authors are partially founded by CONACYT CB-2016-01-284451 grant.

Conflict of interest

The authors declare no conflicts of interest in this paper.

Supplementary

A.   Some inference results

In this Section, we show some results about the inference as the posterior variance of parameters $\theta_0, \theta_1$ and $\theta_3$ in Figure 5 (a)-(c). Since parameters $\theta_2$ and $G_b$ are not practical identifiable for most of the subjects, we omit this information. We show inference results for a possible misclassified IGT patient in Figure 5 (d) and Table 2.

B.   Identifiability

In this section, we use the Similarity Transform Method ^[26] to show some results about the identifiability of parameters $\theta_0$ and $\theta_3$. Let us consider $X = G-G_b$ and the system for $X \geq 0$

subject to glucose observations

We need to find $T$, non singular, such that

After some computations, we can deduce that

then

Since $T$ must be non singular, then $T_{44} \neq 0$ and $\tilde{\theta_0} = \theta_0$. We have $\theta_0$ is identifiable. Now, from the observation condition, we have $T_1 = (1, 0, 0, 0)$, then

Now, $T_{21} = 0$ and $T_{22} = 1$ and $\theta_1$ is identifiable,

Then $T_{23} = 0$ and $T_{33} = 1$

Finally, since $\tilde{\theta_0} = \theta_0$, we have that each

with $T_{44} = \dfrac{\tilde{\theta_3} + \theta_0 (\theta_0 - \lambda_5)}{\theta_3 + \theta_0 (\theta_0 - \lambda_5)}$, is a similarity transform for the system and $\theta_3$ is unidentifiable.

C.   Stability analysis

The process of maintaining blood glucose at a steady-state level is known as glucose homeostasis ^[2]. In this case, we model this steady-state as the basal glucose level, denoted by $G_b$. The model proposed in equations (2.12)-(2.18) has an unique equilibrium point $(G_b, 0, 0, 0, 0, 0, 0)$. This model is piecewise linear, and we would like to demonstrate that every solution of the ODE system converges to this point which in practice means that under a perturbation or stimulus of the glucose level, the body looks for returning to balance. Note that we are not interested in the period of time it takes, even if for diabetic patients we expect greater times than for healthy patients. Let us consider the change of variable $X = G-G_b$. The system given by equations (2.12)-(2.18) may be expressed as

for $x\geq 0$, and by

for $x < 0$. The unique equilibrium point of this piece-wise linear systems is $(0, 0, 0, 0, 0, 0, 0)$. According to ^[35], the switched linear system given by equations (C.1) and (C.2) is Globally Uniform Exponential Stable (GUES) iff it is locally attractive for every switching signal (see Theorem 2.4). Let us recall that the local attractivity property means that all trajectories starting in some neighborhood of the origin converge to the origin. For $x\geq 0$, we can convert system in Equation (C.1) in the third order linear equation

with $f(t) \to 0$. Associated to equation (C.3), we have a third order polynomial

which roots determine the behavior of its solutions. Let us recall that the discriminant of a depressed cubic $t^3 + pt + q = 0$ is $\Delta = -(4p^3 + 27q^2)$. For equation (C.4), after a change of variable to obtain it depressed form, we have

Recall that for $\Delta < 0$, the equation has one real root and two non-real complex conjugate, while for $\Delta > 0$, three real roots. Note that the value $\theta_1 = \dfrac{16}{27}\lambda^2_5$ satisfies $\Delta = 0$.

For $\theta_1 < \dfrac{16}{27}\lambda^2_5$, $\Delta > 0$ is satisfied. In this case is very simple to verify that the three real roots are negative. Consider the polynomial $G(t) = t^3 + 4 \lambda_5 t^2 + 4 \lambda^2_5 t + 2\lambda_5 \theta_1$, which satisfies:

that means, there exist a root $x_1 \in ]-2/3 \lambda_5, 0[$ of $G(t)$. Analogously, we can find $x_2 \in ]-2 \lambda_5, -2/3 \lambda_5[$ and $x_3 \in ]-4 \lambda_5, -2 \lambda_5[$ roots of $G(t)$. Since all are negative, the local attractivity property is satisfied in this case.

For $\theta_1 > \dfrac{16}{27}\lambda^2_5$, $\Delta < 0$ is satisfied. In this case, the roots of $G(t)$ are given by

with $\xi = \dfrac{-1 + \sqrt{3}i}{2}$, $C = \left(\dfrac{\Delta_1 + \sqrt{\Delta_1^2 -4 \Delta_0^3}}{2} \right)^{1/3}$, $\Delta_0 = b^2-3ac$ and $\Delta_1 = 2b^3 - 9abc + 27a^2d$. Since $a = 1, b = 4\lambda_5, c = 4\lambda^2_5, d = 2\lambda_5 \theta_1$, we have

then

Since $\theta_1 > \dfrac{16}{27}\lambda^2_5$ then $C > 0$, and therefore $x_0 = -\dfrac{1}{3a}\left(b + C + \dfrac{\Delta_0}{C}\right)$ is real and negative. Now, $x_1 = -\dfrac{1}{3a}\left(b + \xi C + \dfrac{\Delta_0}{\xi C}\right) = -\dfrac{1}{3a}\left(b + \xi C + \xi^2 \dfrac{\Delta_0}{ C}\right) = -\dfrac{1}{3a}\left(b -\dfrac{1}{2} C - \dfrac{\Delta_0}{2 C} + \dfrac{\sqrt{3}}{2}\ i \left(C - \dfrac{\Delta_0}{ C} \right) \right)$. Its real part is $-\dfrac{1}{3a}\left(b -\dfrac{1}{2} C - \dfrac{\Delta_0}{2 C} \right)$. We are interested in determining under what conditions $\left(b -\dfrac{1}{2} C - \dfrac{\Delta_0}{2 C} \right) > 0$. A straightforward computation leads to the condition

This condition establish an upper bound for the value of $\theta_1$. In practice, this condition is satisfy for values of $\theta_1 < 29$. A similar analysis is valid for $x < 0$ and $\theta_2 < 29$.