Kronecker product decomposition of Boolean matrix with application to topological structure analysis of Boolean networks

1.   Introduction

According to the CDC, nearly 9.4% of the U.S. population has diabetes (30.3 million) including individuals with type 1, type 2 and gestational diabetes ^[1]. Type 2 diabetes mellitus (T2DM) is usually diagnosed in adulthood when the body cannot regulate glucose levels properly without medication or treatment. T2DM accounts for 90–95% of diabetes cases (29.1 million) and manifests in 15–30% of pre-diabetes cases (86 million), a population where 9 out of 10 individuals do not know that they are pre-diabetic ^[2]. In 2017, an estimated 84.1 million adults had prediabetes, a condition defined by higher than average blood glucose levels but not high enough to be diagnosed with diabetes ^[1]. Individuals with prediabetes have a higher risk of developing T2DM and other cardiovascular diseases ^[1]. In fact, an estimated 15–30% of the U.S. population with prediabetes are more likely to develop T2DM within 5 years. This is sobering information that must be understood by a subpopulation where 9 out of 10 individuals are unaware of their pre-diabetic status ^[3,4].

Since T2DM is a chronic disease, the quality of life of patients with diabetes decreases in the long-term depending on severity and other co-existing health conditions. Diabetes increases risk for blindness, kidney failure, high blood pressure, heart disease, stroke, amputations, dental disease, depression, and pregnancy complications ^[2]. Prevalence levels are high among American Indians and Alaska Natives (15.9%), Non-Hispanic blacks (13.2%), and Hispanics (12.8%) when compared to its prevalence among Non-Hispanic whites (7.6%) and Asian Americans (9.0%) ^[5]. Diabetes is unfortunately also correlated with increases in individuals' economic and health burdens, thereby increasing mortality and morbidity risks. The annual costs of diabetes has been estimated to be about 245 billion dollars, 176 billion attributed to direct and 69 billion costs that include (disability, premature death, and work absenteeism) ^[2]. Among diagnosed diabetes cases seeking treatment, 14% use insulin only (2.9 million adults), 56.9% use oral medication only (11.9 million adults), 14.7% use both insulin and oral medication (3.1 million adults), and 14.4% use neither insulin nor oral medication (3 million adults) ^[2].

1.1.   The role of free fatty acids on the progression of diabetes

The clinical manifestation of T2DM is characterized by hyperglycemia (prolonged blood glucose level above normal), hyperinsulinaemia (a condition in which there are excess levels of insulin circulating in the blood) along with insulin resistance. Insulin resistance (IR) is considered an indicator of the inability of cells in muscle, fat, or liver tissue to respond to insulin, which in turn, limits the ability of the individual take up glucose from the blood easily. Hence, the pancreas must produce and secrete additional insulin, which is needed to bring glucose levels back to the normal range. IR can also develop in the liver, known as hepatic insulin resistance, which is defined as impaired insulin signalling affecting the liver. Clinical studies show that insulin resistance in skeletal muscle and adipose tissue are common in individuals with T2DM ^[6,7,8,9], and thus IR is considered a strong predictor for the development of T2DM.

Several studies support the hypothesis that prolonged presence of free fatty acids (FFA) circulating in the bloodstream or in organs (e.g., skeletal muscle, pancreas, or liver) contributes to the development of T2DM by impacting insulin signalling, glucose transport, insulin secretory function, or other mechanisms leading to disruptions in glucose-insulin regulation ^[10,11,12]. FFA is used as fuel for the liver, resting skeletal muscle, renal cortex, and myocardium ^[10]. During starvation or exercise, the demand for FFA as a fuel source increases, where after an overnight fast lipid oxidation can account for over 70% of total body energy expenditure ^[10]. When FFA is needed as a source of energy, then triglycerides (TG) are broken down into FFA and released into the bloodstream via lipolysis. Hence, FFA can be approximated by TG. Since TG is associated with cardiovascular diseases ^[13], then high TG levels coexisting with normal or impaired fasting glucose further increase risk of metabolic-related diseases, such as T2DM. FFA is stored in adipocytes (adipose cells) from TG which are made up of three FFA molecules and one glycerol molecule. As of today, several hypotheses have been proposed characterizing the mechanistic role of FFA on the progression of insulin resistance (summarized in Table 1). It is widely accepted that elevated FFA levels promote insulin resistance in skeletal muscle and liver ^[10,14,15] and yet the mechanisms underlying the progression of T2DM based on obesity levels are not well understood. Here we formulate a mathematical framework based on what is currently known about the role of FFA on IR progression as a starting point, or baseline model, to be expanded upon in future work to test different hypotheses through computational studies.

1.2.   Bariatric surgery as a treatment strategy

Obesity, characterized by excessive body fat, is considered a risk factor for diabetes since among those individuals diagnosed, 84.7% turned out to be overweight, with 56.9% obese ^[1]. In many cases, patients with severe obesity and diagnosed with T2DM can have improved metabolic health following bariatric surgery, which is a treatment strategy for patients with severe obesity. These procedures may include: gastric banding (such as adjustable and nonadjustable bands), gastric bypass (such as Roux-en-Y variations or any other procedure combined with gastric bypass), gastroplasty (such as vertical banded gastriplasty), biliopancreatic diversion or duodenal switch (such as various modifications), or a combination of these methods ^[24]. In a review of 134 studies ^[24], it was found that 76.8% patients who underwent bariatric surgery improved diabetes-related outcomes, and a total of 85.4% patients had either resolved or improved diabetes-related outcomes. Bariatric surgery patients with diabetes have demonstrated recovery of acute insulin response, also referred to as the "first phase" of insulin secretion ^[25], decreases of inflammatory indicators (C-reactive protein and interleukin 6) which are associated with cardiovascular and metabolic health ^[26], improvement in insulin sensitivity ^[27,28], significant changes in the response of enteroglucagon, defined as a peptide hormone that helps to facilitate the absorption of nutrients in the presence of glucose and fats, to glucose levels ^[29], significant reduction in ghrelin, also referred to as the "hunger hormone" which triggers appetite ^[30], and significant improvement in beta-cell function ^[31]. Hence, a reduction of body fat following bariatric surgery leads to improvements in the patient's metabolic health.

In this paper, we consider the case of a prototypic individual with genetic susceptibility to T2DM, that is, an individual who self-identified as having a family member with T2DM or other metabolic-associated diseases, and thus, could be susceptible defined by its inability to regulate FFA or by a disruption on the regulation of FFA. We refer to this process as the "harmful cycle hypothesis" and define the mechanism as follows. A dysfunctional regulation of FFA that leads to an increase demand for insulin, a demand that is stimulated by both glucose and FFA, which in turn, puts the individuals at higher risk of hyperglycemia when beta-cell compensation fails. In the long-term, insulin resistance is observed when insulin does not effectively regulate plasma glucose and FFA. We introduce a model, a hypothesis-driven model, built off of models previously used to study the joint dynamics of insulin, glucose, and FFA. These dynamics are studied quantitatively and qualitatively within a framework that includes model validation. Parameters are fitted to data obtained from a heterogeneous sample of patients ranging from non-diabetic and non-obese to diabetic and obese so that the model results to may provide insights into the physiological factors that govern glucose homeostasis. By fitting the model to clinical data of patients who underwent bariatric surgery, we explore to what extent does a reduction of body fat following surgery alter a patient's metabolic health by recovering insulin sensitivity, beta-cell function, and other related parameters. This paper is organized as follows: In section 2 the mathematical models are described including the classic minimal model (in section 2.1) of insulin, glucose, and FFA as described in the work of ^[32] and an explicit time-delay model of insulin, glucose, and FFA dynamics (in section 2.2) incorporating the harmful cycle hypothesis with corresponding analytical results, in section 4 the results of the model validation and parameter estimates are shown to compare both the minimal model and explicit time delay model, and in section 5 the conclusions and future work is described.

2.   Mathematical modeling framework

The insulin sensitivity is defined through the clamp protocol. Although the clamp protocol is the gold standard for assessing insulin sensitivity, not only do the subjects suffer pain from the procedure, but also the test is labor-intensive and financially expensive. The intravenous glucose tolerance test (IVGTT) has been clinically considered the most accurate protocol next to the clamp protocol to determine insulin sensitivity and glucose effectiveness through the approach of mathematical modeling with curve-fitting ^{[34,35,36,37,38]}. The data provided through an IVGTT offer rich information and offer a more realistic picture of a subjects metabolic portrait of insulin sensitivity and glucose effectiveness by differentiating glucose production from the liver and the exogenous glucose administered during the study protocol ^[39]. Furthermore, the IVGTT protocol has been extended to gain a better understanding of the dynamical regulations underlying FFA, insulin and glucose ^[32,40], in which, as the original IVGTT, subjects fast overnight and then are given a bolus of glucose infusion intravenously (e.g., 0.33 g/kg of body weight or 0.3 g/kg body weight of a 50% solution), which is administered into the antecubital vein in approximately 2 minutes ^[39,41]. FFA, plasma glucose and insulin levels are sampled over the duration of the test.

The short dynamics captured by the mathematical model begins with a rise in plasma glucose due to bolus infusion, which triggers pancreatic beta-cells to quickly secrete insulin into the bloodstream. Insulin mediates glucose removal, also referred to as insulin-dependent removal, which in turn, lowers plasma glucose to basal level and then the demand for insulin is inhibited, i.e. negative feedback. Meanwhile, FFA production is inhibited by insulin when glucose supply is high (see Figure 1). On the other hand, some prior studies suggest that insulin inhibition of FFA is weak in individuals genetically predisposed to metabolic-related diseases. Furthermore, it is known that FFA may enhance basal and glucose-stimulated insulin secretion among individuals including those with diabetes ^[10]. We hypothesize and validate that insulin ineffectively regulates FFA and higher FFA may reduce glucose transport, which leads to a harmful cycle promoting hyperglycemia and developing insulin resistance in the long-term.

2.1.   Minimal model of glucose, insulin, and FFA

The minimal model was introduced in 1979 and 1980 by Bergman, Cobelli and their colleagues. It was the first model to define two significant indices, the glucose effectiveness index and the insulin sensitivity index, which quantify two clinically and physiologically relevant features ^[34,35], and to date, continues to be used in clinical settings or improved in modeling studies. Insulin kinetics, including both first phase and second phase insulin secretion, reflects pancreatic responsiveness. Insight into an individual's glucose tolerance or intolerance are obtained via the estimation of pancreatic responsiveness, glucose disappearance, and insulin sensitivity. The minimal model has been extended and is now widely used in various experimental settings. The software MINMOD based on the minimal model is used by clinicians and researchers who are interested in quantifying insulin sensitivity and beta-cell responsiveness ^[42]. However, fewer mathematical models have been proposed that link insulin, glucose, and FFA ^{[32,40,43,44]}. The model developed by Chow and his colleagues ^[40] captures the interactions between remote insulin $X(t)$, glucose $G(t)$, and FFA $F(t)$. Glucose enters the body intravenously and then it is removed from immediate use by other tissues at a constant glucose effectiveness rate $S_G$ or via insulin-mediated removal modeled by the interaction term $S_I XG$, where $S_I$ represents insulin sensitivity. A proportion of insulin is available for use at rate $c_X$, while the remainder is removed either by natural degradation or by the kidney and liver. The maximal lipolysis rate is given by $l_0 + l_2$. Insulin inhibition of lipolysis is denoted by $X_2$ with the exponent $A$, and the clearance rate of FFA is denoted by $c_f$. The system of equations describing these dynamics is given by

where $I(t)$ represents the insulin concentration in the body over the time $t$. In minimal model and its siblings in ordinary differential equation systems, the known physiological delay of insulin secretion into the body in response to the rise of glucose level is incorporated using compartment-split techniques, where the insulin compartment is split into two equations $I(t)$ and $X'(t)$ leading to higher dimensional systems of ordinary differential equations ^{[34,35,41,45,46,47]}. The Model (2.1) is next modified via the incorporation of an explicit fixed time-delay that is explicitly linked to insulin secretion.

2.2.   Explicit time delay model of glucose, insulin, and FFA

A time delay for the insulin secretion in the body is a key physiological factor in the endocrine regulation of insulin, glucose and FFA, since without insulin is essential for lowering glucose levels and suppressing the production of FFA during the hyperglycemic state ^[7]. This delay can be incorporated explicitly in a system of delay differential equations in order to investigate more realistic intrinsic phenomena in the biological process ^[39,41]. More recent models incorporating an explicit time delay provide more accurate quantification of insulin sensitivity and glucose effectiveness since these models are more robust ^[39,41]. Mathematical analysis and numerical simulations ^[39,48,49] show that these models generate results that can match the observations generated by clinical studies ^[51,50,55], both in the short- and long-term dynamics (see work by ^{[39,48,52,53,54]}). The above model is an extension of the insulin and glucose model studied in ^[39] and the equation representing FFA was adapted from the work in ^[40]. We hypothesize that when insulin ineffectively regulates FFA, then higher FFA may reduce glucose transport, leading to a harmful cycle promoting hyperglycemia and thus contributing to the development of insulin resistance in the long-term.

Next, we consider the interplay of glucose denoted by $G(t)$, insulin denoted by $I(t)$, and FFA denoted by $F(t)$. The parameter $(S_iI_b + S_g)G_b$ represents the rate of constant average hepatic glucose input in the short dynamics, that is, in the beginning of the IVGTT protocol. The elevated glucose is either immediately used from other cells at a rate $S_g$ or by insulin-mediated removal by the interaction term $S_iGI$ at the rate $S_i$ representing insulin sensitivity. Both glucose-stimulated and FFA-stimulated insulin production follow the dynamics of Hill's function, $\frac{\sigma_1G(t-\tau)^\gamma}{\alpha^\gamma + G(t-\tau)^\gamma}$ and $\frac{\sigma_2F(t)^\beta}{\sigma^\beta + F(t)^\beta}$, respectively. Here $\tau$ represents the explicit time delay for glucose-stimulated insulin secretion. Insulin has a natural degradation rate $d_i$. For FFA, $g_0 + g_1$ represents the maximal lipolysis rate. The concentration of insulin inhibition of lipolysis is $I_2$ with the exponent $\kappa$, and the FFA clearance rate is denoted $d_f$. The model is illustrated in Figure 2.

Hence, this system is given by

with positive initial conditions, where the parameters $\beta$ and $\gamma$ represent the Hills function coefficient and $\sigma$ and $\alpha$ represent the values of half-saturation, respectively. The descriptions of the parameters for Model (2.2) are shown in Table 3.

3.   Mathematical treatment of explicit time delay model

Basic properties of the model. In this paper we only show some biologically relevant results of Model (2.2). First we show that the basal levels of glucose, insulin and FFA are the only equilibrium point of the model.

Theorem 3.1. Model (2.2) assumes a unique steady state $E_b = (G_b, I_b, F_b)$, where $G_b, I_b$ and $F_b$ are basal levels of glucose, insulin and FFA, respectively.

Proof. In fact, it is obvious that $(G_b, I_b, F_b)$ is a steady state of Model (2.2). Suppose Model (2.2) has another steady state $(G_*, I_*, F_*)$, a solution of the following system

We will show that $G^* = G_b, I^* = I_b$ and $F^* = F_b$. If $G^* > G_b$, then (3.1) implies that $I^* < I_b$, and thus $F^* > F_b$ by (3.3). Therefore $G^* > G_b$, $F^* > F_b$ and (3.2) lead to a contradiction $I^* > I_b$.

If $G^* < G_b$, then (3.1) implies that $I^* > I_b$, and (3.3) implies that $F^* < F_b$. Again $G^* < G_b$, $F^* < F_b$ and (3.2) lead to a contradiction $I^* < I_b$.

Theorem 3.2. All solutions of Model (2.2) with positive initial conditions are positive and bounded.

Proof. Let $(G(t), I(t), F(t))$ be a solution of Model (2.2) with $G(0) > 0, I(0) > 0$ and $F(0) > 0$. Assume that $G(t)$ is non-positive for some $t$, then there must exist a $t_0 > 0$ such that $G(t_0) = 0$ and $G(t) > 0$ for $0 \leq t < t_0$. Moreover, then $\frac{dG(t_0)}{dt}\leq0$, which contradicts the following

Therefore $G(t) > 0$ for all $t > 0$. Similarly, for $I(t)$, assume that $\exists$ $t_1 > 0$ such that $I(t_1) = 0$ and $I(t) > 0$ for $0\leq t < t_1$. Then, $\frac{dI(t_1)}{dt} \leq 0$, which is a contradiction to

Therefore $I(t)$ is positive for $t > 0$. Finally, the similar and standard approach as above ensures $F(t) > 0$ for $t > 0$.

For the boundedness, clearly $G'(t) \leq (S_iI_b + S_g)G_b - S_g G$, and thus $0 \leq G \leq {(S_iI_b + S_g)G_b}/{S_g}$.

Similarly, it can be shown that $F \leq \frac{g_0 + g_1}{d_f }$, and $I(t) < (\sigma_1+\sigma_2)/d_i$ for $t > 0$ by standard treatments.

Remark 1. Theorems 1 and 2 assure that the system of equations in Model (2.2) is well posed, that is, it supports positive bounded solutions under any positive initial condition.

3.1.   Equilibria and local stability

The steady state is obtained by setting equations in Model 2.2 equal to 0 shown in Eqs (3.1)–(3.3). Rearranging terms and substituting $G^*$ and $F^*$ into Eq (3.2), yields the equilibrium point implicitly in terms of $I^*$:

In order to determine the number of roots in the system, let us consider:

It is easy to see that $y_1(I^*)$ and $y_2(I^*)$ must intersect once, which indicates that (3.6) has at least one positive root. Further, we can show that $y_1'(I*) < 0$, and therefore the system (2.2) has a unique steady state. Let us define $z_1(I^*) = \left(\frac{(S_iI_b + S_g)G_b}{S_iI^* + S_g} \right)^\gamma$ and $z_2(I^*) = \left(\frac{1}{d_f } \left(g_0 + \frac{g_1}{1+\left(\frac{I^*}{I_2}\right)^\kappa} \right) \right)^\beta$. Substituting $z_1(I^*)$ and $z_2(I^*)$ into (3.6) yields,

Next, obtaining $\frac{d}{dI^*}y_1(I^*)$ gives,

where

and

It is easy to see that $\frac{d}{dI^*}y_1(I^*) < 0$, and since $\sigma_2 < < \sigma_1$, then this statement is true and there is one unique steady state. This conclusion is consistent with findings of other glucose-insulin regulation models in the literature and demonstrate that our model is well-posed ^{[32,34,35,39,40,41,43,44]}.

We now turn to study the stability of the unique steady state. It is straightforward that the characteristic equation of Model (2.2) is given by

where $b_1 = S_i I_b + S_g, b_2 = d_i+d_f, b_3 = d_id_f - \hat{B}\hat{C}, b_4 = S_i G_b \hat{A},$ $\hat{A} = -\frac{\alpha^\gamma \gamma \sigma_1 G_b^{\gamma -1}}{(\alpha^\gamma + G_b^{\gamma})^2}, \hat{B} = -\frac{\sigma_2 \beta F_b^{\beta - 1}\sigma^\beta}{(\sigma^\beta + F_b^\beta)^2}, ~~ and ~~ \hat{C} = \frac{g_1 \kappa \left(\frac{I_b}{I_2}\right)^{\kappa -1}}{I_2\left(1 + \left(\frac{I_b}{I_2}\right)^\kappa \right)^2}.$

Case with no delay. We analyze the local stability of a positive equilibrium point $E_b$ for our system with no time delay by evaluating Eq (3.7) with $\tau = 0$, which gives

Applying Routh-Hurwitz Stability Criterion ^[59] for a cubic polynomial, that is, for the cubic polynomial:

where all $a_i$ are positive. The Routh array is

so the condition that all roots have negative real parts is $a_1a_2 > a_0a_3$. Therefore, in this case, the equilibrium point is asymptotically stable if

Case with delay. Next, we investigate the stability for the case $\tau > 0$. Notice that the equilibrium point $E_b$ is stable for $\tau = 0$. If there is some $\tau > 0$ such that $E_b$ is unstable, then the characteristic Eq (3.7) must have a pair of pure imaginary roots $\pm\omega i$ with $\omega > 0$ ^{[48,49,56,57]}. Thus the characteristic equation becomes

After algebraic rearrangement,

Then we have,

which leads to

Let $u = \omega^2 > 0$, then

where $c_2 = [(b_1+b_2)^2-2(b_1 b_2 + b_3)]$, $c_1 = [(b_1 b_2 + b_3)^2-2b_1 b_3 (b_1 + b_2)-b^2_4]$, and $c_0 = [(b_1b_3)^2-(b_4d_f)^2]$.

According to Descartes' Rules, Eq (3.10) has no positive root if $c_2, c_1, c_0 > 0$. Numerically, we test these conditions and observe that the conditions are met for all groups (see Table 4). Therefore we summarize the above results in the following.

Theorem 3.3. If $c_2, c_1, c_0 > 0$, the unique equilibrium point $E_b$ is asymptotically stable.

4.   Results of model validation and parameter estimation

To determine whether Model (2.2) could capture the metabolic profile of an individual, we fit the model to the clinical data in order to estimate the physiologically relevant parameters. Data were obtained from a study ^[58], where an IVGTT was given to N = 38 patients (14 men and 24 women) with severe obesity who underwent bariatric surgery and varied metabolically. The control group was neither obese nor diabetic and did not undergo bariatric surgery. It consisted of 6 men and 6 women (N = 12). The treatment group consisted of 14 men and 24 women (N = 38) and included severely obese individuals divided into 3 groups: First, the normal fasting glucose (NFG) tolerance group (N = 9); second, the impaired fasting glucose (IFG) group (N = 17); and third, the T2DM group (N = 12). Participants in the control and NFG groups were in similar age groups: 36.7 $\pm$ 1.9 and 35.9 $\pm$ 3.4 years, respectively. Similarly, the IFG and T2DM groups were closer in age: 45.2 $\pm$ 2.5 and 44.6 $\pm$ 2.4, respectively. The body mass index (BMI) calculated as weight divided by height squared ($kg/m^2$) were significantly reduced post-surgery. The control group had an average BMI of 23.1 $\pm$ 0.7 at baseline. Participants in the NFG, IFG, and T2DM groups underwent bariatric surgery and were considered to have severe obesity (BMI $> 40 kg/m^2$). The average BMI decreased from baseline levels in the: NFG group from 48.6 $\pm$ 1.7 to 32.4 $\pm$ 1.4, IFG group from 58.1 $\pm$ 1.4 to 39.6 $\pm$ 1.4, and T2DM group from 53.9 $\pm$ 1.7 to 39.6 $\pm$ 1.5 after surgery.

Plasma glucose, insulin, and FFA levels were measured by the IVGTT before and seven months after undergoing bariatric surgery (N = 38). The control group did not undergo bariatric surgery and were sampled at baseline and at follow-up in seven months. The IVGTT study protocol was performed after a 10 to 12 hours fast. Baseline blood samples were obtained before glucose administration at 15, 10, and 5 min time marks. At time 0 glucose was administered (50% dextrose; 11.4 g/m2 body surface area) in less than one min. Blood samples were obtained after 2, 3, 4, 5, 6, 8, 10, 12, 14, 16, 19, 22, 25, 30, 40, 50, 60, 70, 80, 100, 120, 140, 160, and 180 min to measure glucose, insulin, and FFA concentration levels ^[58]. Here the minimal model (Model (2.1)) and the explicit time delay model (Model (2.2)) are fit to a dataset ^[58] and model parameters were estimated using the function fmincon for nonlinear programming in MATLAB R2019a. The data are the average values of multiple individual data. The data was extracted from the paper using "Plot Digitizer." Initial conditions from the data are considered for time $t = 0$.

4.1.   Validation of the explicit time delay model

The performance of each model on predicting the qualitative trends of the clinical data were determined through many computational studies and numerical simulations. The estimated model parameters for the explicit time delay model is shown in Table 5 and simulations are shown in Figure 3. Qualitatively, comparing the plasma glucose levels pre- and post-bariatric surgery, show a significant improvement, where the NFG, IFG, and T2DM data nearly resembles the glucose trends of the control group. The overall glucose trends for each group pre- and post-surgery were captured by Model (2.2) (shown in Figure 3). The constant rate of glucose effectiveness ($S_g$) and constant insulin sensitivity rate ($S_i$) had insignificant changes in the control group but significantly improved (e.g. increased) for the NFG, IFG, and T2DM groups.

The plasma insulin levels significantly improved comparing pre- and post-bariatric surgery results, where the NFG, IFG, and T2DM data closely match the overall trends of the control group. Model (2.2) captures the overall trends for insulin (shown in Figure 3). The constant insulin degradation rate ($d_i$) decreased from post- compared to pre-surgery for all groups and remained constant for the control group. The maximum secretion rate ($\sigma_1$) decreased post-surgery for IFG and T2DM but not for NFG nor control groups. The constant rate of FFA-stimulated insulin secretion ($\sigma_2$) increased for all groups except for T2DM.

We also evaluated changes in the estimated explicit time delay parameter $\tau$. As expected, the time delay in the control group pre- and post-surgery were similar since they did not have surgery. The estimated time delay decreased after surgery in the NFG and IFG groups, which supports prior findings that have observed improvements in acute insulin response and beta-cell function after bariatric surgery ^[25,31]. In contrast, an overall increased time delay was observed in patients with T2DM. One possible explanation is that the treatment strategy implemented for changes and/or the presence of other clinical health factors may explain variations in these metabolic parameters.

The variance for plasma FFA levels significantly reduced comparing pre- and post-bariatric surgery results, where the NFG, IFG, and T2DM data nearly overlap. The model results qualitatively match the overall trends, except for the NFG group (shown in Figure 3). The maximal lipolysis rate ($g_0 + g_1$) reduced significantly from pre- to post-surgery for the IFG and T2DM groups. The insulin inhibition ($I_2$) rate also decreased post-surgery compared to pre-surgery in the NFG, IFG, and T2DM groups, whereas no changes in FFA clearance rate ($d_f$) were observed.

4.2.   Comparison of minimal model and explicit time delay model

Model validation and parameter estimation were completed for Model (2.1) and Model (2.2) with the estimated model parameters for each model shown in Tables 5 and 6, respectively. The numerical results are shown in Figure 3. Both models qualitatively predicted the overall trends for the control, NFG, IFG, and T2DM groups of plasma glucose. In this case, both models fit the data well. Overall, both Model (2.1) and (2.2) fit the insulin data adequately well for all groups. However, the steady state values produced by Model (2.1) underestimated solutions in comparison to the actual data for the control, NFG, IFG, and T2DM groups both pre- and post-surgery. In the control, NFG, and T2DM groups, the solutions of Model (2.1) falls below the actual data; and the results of Model (2.2) matches the data much better. For the IFG group, Model (2.2) performs much better with describing the insulin dynamics pre- and post-surgery, whereas Model (2.1) does not closely match the data. Qualitatively, the trends were not closely matched by Model (2.1) pre-surgery for the NFG, IFG, and T2DM groups. In this case, Model (2.2) approximated the overall trends more accurately than Model (2.1).

Model validation for plasma FFA levels indicate that the overall trends for FFA were captured slightly better using Model (2.2) for the Control, NFG, IFG and T2DM groups. It also appears that, qualitatively, the FFA dynamics vary more in the transient phase for Model (2.1) compared to Model (2.2). For the control group, Model (2.1) reached the true steady states better compared to Model (2.2). In the NFG group, Model (2.1) fitted the overall trends better than the Model (2.2). However, for the IFG group, Model (2.2) fitted the data better than Model (2.1). Similarly, in the T2DM group, Model (2.2) also fitted the data better than Model (2.1).

To determine the overall goodness of fit, the Akaike information criterion (AIC) was calculated (see Table 7). The AIC allows us to determine the overall quality of the statistical model for a given set of data and estimated parameter values. The delay model performed better for the NFG (Total AIC = 1014.47), IFG (Total AIC = 811.09), and T2D (Total AIC = 799.56) group settings in comparison to the minimal model (see Table 7). In summary, the AIC for all groups combined was lower for the delay model (AIC = 2625.12) than in the minimal model (AIC = 2910.85).

5.   Discussion

We investigated the efficacy of insulin suppression on lipolysis and assessed the hypothesis that FFA-stimulated insulin secretion might play a vital role on the progression of insulin resistance. While the role of FFA-stimulated insulin secretion on the progression of insulin resistance is known, we formulated a mathematical model that incorporates FFA-stimulated insulin secretion explicitly and were able to quantify differences in this effect pre- and post-surgery. In Table 5, we found that $\sigma_1$ and $\beta$ estimated pre-surgery is much greater than the estimated values post-surgery in the impaired fasting glucose (IFG) and T2DM settings. Hence, our findings indicate that the effect of FFA-stimulated insulin secretion on insulin and glucose levels in patients with IFG and T2DM resembled the estimates obtained for the control group indicating that the bariatric surgery patients could recover moderate insulin secretion activity stimulated by FFA. Similarly, the efficacy of insulin inhibition of FFA. Similarly, $I_2$, which represents the efficacy of insulin inhibition of FFA production, on the regulation of FFA, insulin and glucose in patients with NFG, IFG, or T2DM resembled the control group after surgery. These findings indicate that the reduction in weight or body mass due to surgery improve insulin action. This reinforces what we know, that insulin action is essential for regulating plasma FFA and glucose levels. Our results highlight the robust effect of insulin, that is, it can be observed not only in the short-term dynamics during an IVGTT protocol but also in the long-term in longitudinal epidemiological studies.

Overall, it seems that the inclusion of the explicit time delay in Model (2.2) was able to capture the overall qualitative trends well except for the NFG group when compared to Model (2.1). Particularly, Models (2.1) and (2.2) captured the peak and steady states for glucose concentration levels pre- and post-surgery in all groups. Both models captured the initial peak in insulin levels as well as the steady states for the control group. Though the steady states and rapid decline of insulin levels predicted by the models were not exactly equivalent to the averaged data for the NFG group, both models were able to predict the rise and decline in insulin with (Model (2.1)) and without (Model (2.2)) time delay. As shown for the IFG group, the explicit time delay in Model (2.2) captured the two peaks in insulin, which has been shown in previous studies ^[39,41]. Further, both models captured overall rise, decline, and steady state for insulin levels in the IFG and T2D groups. In contrast, the models fit FFA levels over time less precisely. Among the control group, the initial high basal levels and decrease in FFA were predicted by both models. Among the IFG, NFG, and T2D groups the fit to the FFA data were less precise. However, both models had an initial decline and rise of FFA levels until reaching steady states. Overall, obtaining the final steady states for FFA were difficult to compute across all groups due to the complexity of the model as well as large variance in study participants that could be driven by several biological factors, such as age, gender, stage of diabetes or diabetes-related conditions, and other factors.

The results show that the incorporation of an explicit time delay in the model led to better approximations of the qualitative dynamics when compared to those generated by the minimal model. In previous studies ^[22,38], a key feature of $\tau$ is not only the biological interpretation given by quantitatively estimating the fixed time delays across patients and comparing their differences, but also its ability to capture different qualitative behaviors (or peaks) that may give insight into beta-cell function among patients with and without impaired fasting glucose. In this study, it was shown that the explicit time delay model could capture the two peaks for the averaged data of the impaired fasting glucose group, but the interpretation of $\tau$ for a sample average of this group and between groups is limited due to the data being averaged since the proposed model is intended to model the dynamics of a single person, and not a group of individuals. In our future work, we plan to analyze individual patient data and explore the heterogeneity in estimates of $\tau$ across different groups based on metabolic health, gender, and age.

The novelty of our work is proposing a dynamic model of glucose, insulin, and FFA with the inclusion of an explicit fixed time delay. Additionally, our model was constructed to directly study the harmful cycle that is activated when FFA is not regulated properly by insulin and, in turn, FFA continues to circulate in the bloodstream, eventually promoting insulin resistance. Though we used averaged data for our model validation, our computational studies demonstrated that the model qualitatively captured the glucose, insulin, and FFA levels trends over time after an IVGTT protocol. However, the use of average data instead of individual patient data limits the conclusions of this model, and hindered our ability to measure insulin sensitivity for each individual, which is needed to understand and assess insulin suppression by FFA. Individual patient data is needed to better assess models' results since the use of averages eliminates the variance that needs to be quantified in each group. Further work is needed to fit more data to evaluate this observation in depth at the patient level. That is, the variabilities in glucose effectiveness, insulin suppression, and insulin sensitivity cannot be assessed at the individual level using average data between groups, and therefore, patient level data is needed. More specifically, the dynamics of insulin, glucose, and FFA need to be better understood for individuals within and across metabolic groups (e.g., nondiabetic or prediabetic, type 2 diabetic) in order to understand the conditions in which FFA is suppressed effectively and in which cases it is not, and how bariatric surgery impacts the FFA dynamics. Lastly, the physiological parameters adjusted to fit the data matched findings from the literature on bariatric surgery, an indicator that the model captures some of the observed phenomena. Future work would require individual patient data and possibly a new model that better captures the impact of bariatric surgery on insulin regulation of glucose and FFA. Furthermore, while the results of our analyses were similar to the work shown in ^[40], our work provided new insight by the inclusion of $\tau$ in our model as well as other terms to explore the effect of FFA on insulin action. Additionally, models modified to include specific mechanisms corresponding to different bariatric surgery types would also provide greater insight into the roles of FFA and how it is impacted by bariatric surgery.

Acknowledgments

The authors would like to acknowledge Dr. Sherry Towers for her help and suggestions on the technical and statistical aspects of our final manuscript. This project has been partially supported by grants from the Alfred P. Sloan Foundation, the Graduate Assistance in Areas of National Need (GAANN), the National Science Foundation (DMS 1263374), the National Security Agency (H98230- 15-1-0021), the Office of the President of ASU, and the Office of the Provost at ASU. ALM would like to acknowledge funding by the National Heart, Lung, and Blood Institute (NHLBI) of the National Institutes of Health (NIH) grant number T32HL072757. This research is partially supported by the Brown University Data Science Initiative.

Conflict of interest

No potential conflict of interest was reported by the authors.