Glucocorticoids (GCs) are a class of steroid hormones that have profound effects on energy mobilization, especially glucose metabolism. Synthetic GCs are widely prescribed in medical practice because of their anti-inflammatory, immunosuppressive, and antiallergic effects. However, excess and/or long-term treatment of GCs can induce undesired diabetogenic side effects (21). Aside from drug-induced hypercortisolism, pituitary tumor (Cushing’s disease) and chronic stress are the other two causes of excessive GCs and increase the risk of diabetes development (8). Investigations of the links between glucocorticoid and glucose dynamics are desired to achieve effective glucose control.

Glucocorticoids facilitate the process of gluconeogenesis in the liver, while they reduce glucose uptake and utilization by antagonizing insulin effects in white adipose tissue and skeletal muscle. As a result, over-exposure to GCs leads to hyperglycemia and insulin resistance (22). Although the causal relationship between GCs and dysglycemia is affirmative, the impact of GCs on the pancreatic beta-cells remains debatable (21). Several studies proposed that the effects of synthetic GCs on pancreatic islets and insulin biosynthesis or release depend on the dose and duration of GCs treatment (23, 24). Research in murine models and human studies have shown excess GCs can cause compensatory beta-cell hyperplasia and hyperinsulinemia, with the coexistence of normoglycemia. However, long-term GCs therapy that oversteps the beta-cell compensatory capacity begets impaired insulin secretion, hyperglycemia, and consequent type 2 diabetes (22).

Zavala et al. developed a mathematical model investigating the impact of disrupted cortisol rhythms on the response to oral glucose tolerance tests (OGTT) (25). The model incorporates the effect of transmembrane glucose transporters (GLUTs) on the glucose uptake of fat and skeletal muscle cells, which is under the regulation of both insulin and GCs. In particular, insulin facilitates the translocation of GLUT1, GLUT3 and GLUT4 from intracellular pools to the cell membrane to amplify the glucose uptake in adipocytes and muscle cells, while GCs antagonize this process by translocating GLUTs from the cell membrane back to intracellular compartments (26, 27). In pancreatic beta cells, GLUT1, GLUT2 and GLUT3 are involved in glucose sensing and possess different affinities for extracellular glucose compared to those in fat and muscle cells (26, 28). The model is described by the following ordinary differential equations:

Where G(t) and I(t) stand for the blood concentrations of glucose and insulin at time t (min), respectively; the variable T∈(0,1) denotes the fraction of translocatable GLUTs in the cell membrane of peripheral cells, and 1–T represents the fraction of GLUTs that docked inside the cell. The functions in Eq. 1-3 are supported by sigmoidal functions with the general form ϕ ( x , k m , h ) = x h x h + k m h , where k_m stands for the half maximum constant and h is the Hill coefficient. In the glucose equation, F(t) represents the glucose boluses from feeding or OGTTs. The term v_fe (G) stands for the endogenous glucose production from gluconeogenesis and glycogenolysis, and v denotes the maximum rate of the process. The term at [c_Lf_L (G)+ c_Mf_M (G)]T represents the glucose uptake by fat and muscle cells, which depends on the fraction of active GLUTs. The factors c_Lf_L (G) and c_Mf_M (G) stand for the glucose transport mediated by GLUT 1,3 and GLUT 4 respectively, which have different affinities for extracellular glucose. The last term r_GG represents the first-order glucose removal. In theinsulin equation, ε stands for the basal insulin secretion rate; S_β (G) denotes the glucose sensing in beta-cells; σ represents the maximum insulin secretory rate; h_Q (G) stands for the regulatory effects of GCs on beta-cell insulin secretion; r_II represents the first order insulin removal. In the last equation, the term (u+v_If_I (I))(1–T) accounts for the translocation rate of GLUTs to the cell membrane, where u denotes the basal translocation rate and f_I (I) represents the insulin-mediated translocation from the intracellular pools to the cell membrane, at a maximum rate v ₁. In the second term, (d+ v_Qf_Q (Q))T accounts for the translocation rate of GLUTs from the cell membrane down to intracellular pools, where d denotes the basal translocation rate and v_Qf_Q (Q) represents the GCs regulated translocation in the same direction.

The model predicts a magnified glucose and insulin non-oscillatory OGTT response under the sub-chronic treatment of dexamethasone (cortisol agonist). In comparison, excess cortisol may enhance the magnitude of the glucose responses to OGTT with maintained circadian and ultradian variability, while greatly suppressing the insulin response and its circadian and ultradian variability. The results also show that excess cortisol results in a right shift of the Starling’s curve toward higher fasting glucose levels, which reveals the impact of hypercortisolism on the progression of diabetes. Overall, this work illustrates how mathematical modeling can provide circadian timing approaches to interpret clinical data, and the potential of mathematical modeling to facilitate the design of the chronotherapies for diabetes secondary to hypercortisolism.

Epinephrine, also known as adrenaline, is a stress hormone that can cause an acceleration in heart rate and glucose metabolism, as well as an increase in blood pressure and muscle strength. The release of epinephrine is typically increased under acute stress to prepare the body for fight-or-flight response (29). As the link between stress and diabetes progression has been brought to the forefront, the impact of epinephrine on glucose metabolism arouses increasing attention from researchers ( Figure 2 ). In particular, studies have shown that surgery-induced metabolic stress, which causes the acute elevation of epinephrine, can significantly increase the rate of dysglycemia and the mortality of hospitalized patients (30). In type 1 diabetes, the counterregulatory response of epinephrine to hypoglycemia is an important factor to be considered for the design of glucose control strategies (34, 35).

Several mathematical models have been developed to investigate the quantitative influence of epinephrine on the glucose regulatory system. However, the biological mechanisms some of these models built upon are challenged by further findings of the different metabolic effects of short-term versus long-term epinephrine. Many studies list epinephrine as raising blood glucose, inducing potent insulin-counteracting effects when administered in a short term. Nevertheless, chronic infusion of low-dose epinephrine can enhance insulin sensitivity and glucose uptake in skeletal muscle. The divergent impact may be attributable to the acute versus chronic effects of epinephrine on β ₂ stimulation (31). Chronic epinephrine infusion enhanced glycogen synthesis activation and insulin-dependent glucose uptake in rat skeletal muscles (33). Moreover, experiments have shown low dose of β ₂-adrenoceptor agonists can improve glucose tolerance in diet-induced obese mice within 4 days of treatment. Prolonged treatment with the low dose of β ₂-adrenoceptor agonists can further enhance whole-body insulin sensitivity, immensely reduce hepatic glycogen levels, and lower blood glucose levels (32). Assuming the effects of epinephrine on glucose regulation in animal studies are consistent with human dynamics, a mathematical model should be carefully designed in view of the dosage and duration of epinephrine. We propose that models investigating the impact of epinephrine on the short-term glucose-insulin dynamics should be formulated upon different mechanisms, compared to the models studying the effect on the long-term glucose regulation.

Mohammed et al. established a model with the variables of glucose, insulin, beta-cell mass, and epinephrine to study the glucose regulation under the influence of trauma, excitement and/or stress (36). This model, as shown below, was formulated upon the model of Topp et al. (37), adding the variable of epinephrine to the glucose and insulin equations:

where G (mg/dl), I (μU/ml), β (mg) represent the blood glucose concentration, insulin concentration, and the mass of functional beta-cells at time t (days), respectively. The parameter R ₀ stands for the net rate of glucose production per day. The term S_GOG represents insulin-independent uptake of glucose, while S_IIG depicts the insulin-dependent uptake of glucose. In particular, the coefficient S_I (ml/μU/day) represents insulin sensitivity. The insulin secretion from beta-cells is hypothesized to be stimulated by elevated glucose levels in the form of the Hill function, and the parameter σ denotes the secretory capacity per beta-cell. The parameter k is the insulin clearance rate (/day). The functional beta-cell mass is designed as a second degree polynomial function of glucose with the assumption that moderate glucose level facilitates the growth of beta-cells, while high glucose level aggravates beta-cell apoptosis. The term G_e (mg/dl/day) stands for increasing rate of glucose concentration due to epinephrine, and ρ (/day) represents the rate constant of insulin suppression by epinephrine.

As we can deduce from the formulation of the model, Mohammed et al. assumed that epinephrine can raise glucose levels by increasing hepatic glucose production and suppressing insulin secretion. Considering this model was built on Topp’s model, which was developed to study the long-term glucose-insulin dynamics, we expect a better formulation of the model should take the long-term effect of epinephrine into consideration. The enhanced glucose uptake in skeletal muscle and improved insulin sensitivity underthe durable impact of epinephrine can be revealed by adjusting the settings of the model equations. Furthermore, as the long period effect of low-dose epinephrine may reduce the blood glucose level, it seems improper to investigate trauma/stress-induced dysglycemia, focusing only on the impact of epinephrine. For example, exercise, as one of the most potent stimuli to release epinephrine, is a significant approach to ameliorate/prevent diabetes (31). Durable excessive secretion of GCs induced by chronic stress may be a more reasonable factor accounting for the commencement of dysglycemia.

Kwach et al. proposed a model (38) to study the acute influence of epinephrine on short-term glucose-insulin dynamics. The model was built upon the assumption that epinephrine can induce small net stimulation of insulin secretion from pancreatic β-cells, which remains debatable as human studies have confirmed the effect of epinephrine on repressing endogenous insulin secretion (39, 40). The assumed self-stimulating effect of epinephrine, presented in the epinephrine equation of the model, needs further justification as well. Kumar and Sandhya filled in biological details of the model of Kwach et al. in their work (41). Their paper cited the study of Sherwin et al. (42), in which the experiment showed that a rapid riseof epinephrine can induce a transient elevation of hepatic glucose output, suppress endogenous insulin secretion, and directly inhibit insulin-stimulated glucose utilization. Nevertheless, the model of Kwach et al. was directly employed in (41) without further modification. The sign of the epinephrine term might be negative in the insulin equation of the model to describe the negative impact of epinephrine on insulin secretion. Examining the biological mechanisms underlying mathematical models isa crucial step to obtaining constructive model implications. A well-developed model for the acute influence of epinephrine on glucose regulation may help to devise a glucose-control strategy for hospitalized patients at risk of hyperglycemia due to surgery-induced metabolic stress.

Type 1 diabetic patients with exogenous insulin therapy are exposed to the risk of hypoglycemia, as their systemic insulin levels may not be reduced in time when the glucose levels begin to decline. In this case, epinephrine becomes the first line of counterregulatory hormone responding to hypoglycemia, due to the early deterioration of glucagon secretion. Moscardó et al. formulated a model to investigate the counterregulatory action of epinephrine during hypoglycemia in type 1 diabetes (43). The model is extended upon the Bergman Minimal Model (44) by adding terms to the glucose equation to account for the influence of epinephrine:

Where A(0) = A_basal , G(0) = G ^*. The variables G, I, and X stand for the blood glucose concentration (mg/dl), plasma insulin (μU/ml), and insulin action from a remote compartment (/min), respectively. As the insulin supply for type 1 diabetic patients comes from exogenous infusion, the insulin secretion term is denoted by u 1 ( t ) V o l I , where u ₁(t) (µU/min) denotes the insulin infusion rate and Vol_I denotes the insulin distribution volume. The parameter p ₃ represents the rate of insulin input in the compartment. The clearance rates of plasma insulin and insulin from a remote compartment are denoted by n and p ₂ respectively. In the glucose equation, p ₄ stands for the hepatic glucose production rate; u 2 ( t ) V o l G depicts the clamp experimental conditions with the glucose infusion rate u ₂(t) and the glucose distribution volumes Vol_G ; p ₁ G denotes the glucose uptake independent of insulin and XG stands for the glucose uptake rate under the influence of insulin; A _m(t) represents the plasma epinephrine concentration (ng/l) and A_basal denotes the basal epinephrine concentration; G_th represents the glucose threshold activating the epinephrine response; thus, A(t) denotes the epinephrine “effect” which would be activated only after the glucose falls below the glucose threshold; p_hmax(G_{b 2}–G,0) denotes the increase of glucose utilization when glucose stays below a threshold G_{b 2}, the value of which is assumed to be in the hypoglycaemic range and greater than G_th . The term p_hmax(G_{b 2}–G,0) is added mainly for a good fit to data.

Studies have shown some diabetic patients under high insulin therapy experienced the stage where the counterregulatory response of epinephrine begins to prevent plasma glucose from further decreasing after the glucose level falls below the activation threshold of epinephrine (around 60 mg/dl) (43). After the epinephrine concentration peaks at the hypoglycaemic plateau, epinephrine secretion rapidly declines when the glucose level starts to recover, returning to its basal level. The work of Moscardó et al. made it possible to present the physiological behavior during hypoglycemia. As this model neglected the inhibitory effect of epinephrine on insulin-dependent glucose utilization in the short term, an improvement over this model may depict a better influence of epinephrine on short-term glucose regulation. Overall, this modeling approach provides a better understanding of the counterregulatory response of epinephrine and may facilitate the design of predictive methods to avoid hypoglycaemic events.

Growth hormone (GH) or somatotropin, is an important peptide hormone that stimulates growth, cell reproduction, and regeneration. GH also stimulates the production of insulin-like growth factor 1 (IGF-1), a hormone similar to insulin in molecular structure (45). GH has been characterized as one of the anti-insulin hormones. The diabetogenic effect of GH is also validated by the high prevalence of diabetes in patients with acromegaly, a condition wherein the excessive growth hormones are produced citep (29). The risk of hyperglycemia exists likewise in patients with GH deficiency who need the GH administration for treatment. Large-scale cohort studies have shown that, compared to the general population, the incidence of developing type 2 diabetes for children under GH treatment was increased more than six times, especially in patients with predisposing risk factors, such as obesity (46–48).

The effects of GH on glucose regulation ( Figure 3 ) are intricate partially due to its indirect impacts via IGF-1, which has glucose-lowering functions analogous to insulin (49). GH can elevate glucose production in the skeletal muscle and liver and diminish glucose utilization in adipose tissue by antagonizing the action of insulin. Insulin secretion is also enhanced for the compensation of elevated blood glucose after GH administration (50). Prior studies have demonstrated that low-dose GH may have beneficial effects on insulin resistance and glucose homeostasis due to increased circulating IGF-1, while long-term GH treatment in high doses impairs insulin sensitivity and exacerbates insulin resistance (51). Thus, cautious monitoring of the possible adverse impact on glucose metabolism induced by GH treatment is advocated. Mathematical models investigating the quantitative influence of growth hormone on the long-term glucose dynamics may facilitate the examination of the durable effect of GH therapy.

Alali et al. developed a model studying the effect of growth hormone on glucose homeostasis (52). The model, as shown below, is extended upon the model of Boutayeb et al. (53), including an equation of GH to depict the interaction of GH withglucose and free fatty acids (FFA):

Where G, I, β, R, and F stand for the blood glucose level (g/l), plasma insulin concentration (μU/ml), β-cell mass (mg), the fraction of insulin receptors on the membrane of the muscle cells, and the concentration of FFA (μmol/l) at time t (days), respectively. The parameter a denotes the constant secretion rate of glucose by the liver and kidneys. The term b + cRI denotes the total body glucose utilization rate. The FFA-induced glucose production is represented by m ₁(F–F_b ), and the GH-induced glucose production through gluconeogenesis and glycogenolysis is denoted by cGH. In the insulin equation, the factor G 2 ( e + G 2 ) depicts the sigmoidal relationship between the extracellular glucose concentration and the insulin secretion; d 1 + R represents the insulin secretion factor per β cell; fI stands for the insulin clearance by liver and kidneys; fRI denotes the insulin clearance by the muscle cell receptors. The β-cell equation follows the same formulation in Topp’s model (37). In the receptor equation, j represents the recycling rate of internalized receptors: k stands for the insulin-induced down-regulation rate of receptors on the cell membrane; l denotes the clearance rate of the surface receptors. The parameters F_b , G_b , and GH_b represent the basal concentration of FFA, glucose, and growth hormones, respectively. The term m ₃(G–G_b ) denotes the lipogenesis rate induced by excess glucose, and x(GH–GH_b ) represents the lipolysis rate stimulated by GH. The clearance rate of FFA is denoted by m ₂(F–F_b ). In the GH equation, p represents the production rate of GH by the somatotropic cells; ωGH stands for the clearance of GH by the liver; the GH uptake by fat cells and receptors is represented by s(F–F_b ) and zR, respectively.

This work provides the first mathematical model that incorporates GH into the glucose regulatory system. The receptor equation obtained from the base model in (53) impedes this model from data fitting and applications. In addition, the insulin equation cited from the base model may be questionable. The variable R was designed to denote the fraction of insulin receptors on the membrane of the muscle cells, while the term d 1 + R was claimed to represent the insulin secretion factor per eta -cell. The impact of the changed behavior of insulin receptor can be implied by designing the insulin sensitivity coefficient c to vary with other variables in the system. Moreover, as β-cell function and insulin level play a significant role in the glucose regulation, incorporating the variable GH into the insulin and β-cell equation is important for a better investigation of the influence of GH on glucose dynamics. Further modeling work studying the interaction between GH and glucose regulation is desired to help unravel the intricate physiological effect of GH with different dosages and treatment duration.

The hyperthyroid-diabetes model without hyperthyroid treatment is given by

where a ₁(t) stands for the time-dependent FT4 synthesis factor; b ₁ represents the TSH concentration corresponding to half the maximal synthesis rate of FT4; d ₁ is the decay rate of FT4; a ₂ represents the default release rate of TSH from the pituitary when FT4 reaches the euthyroid set point value; b ₂+2U is the concentration of FT4 resulting in half the maximal inhibitory effect on the TSH secretion rate controlled by the pituitary; d ₂ represents the decay rate of TSH. The expression of a ₁(t) varies among individuals with different progression rates of hyperthyroidism.

By setting the coefficients of X in p_i (X) to be close to zero, our hyperthyroid diabetes model can illustrate the case where diabetic patients are free from hyperthyroidism. In the clinical scenario, some patients with hyperthyroidism develop diabetes over time, while others can stay away from diabetes in life time. The diverse genetic traits of individuals may determine their cellular response to the hyperthyroid factor, which can be expressed by different parameter values in p_i (X).

It has long been recognized that the treatment of hyperthyroidism can improve glucose control. A nation-wide cohort study shows that the treatment of thyroid dysfunctions can reduce the manifestation of T2DM (69). To quantitatively analyze the benefits of hyperthyroid treatment to the glucose regulatory system, we incorporate the drug treatment to the FT4 equation and formulate the hyperthyroid-diabetes model under treatment as follows:

where D represents the anti-thyroid drug (ATD, e.g. Carbimazole) dosage (mg); D ₃ stands for the maximum reduction rate of FT4 caused by drug intake; IC ₅₀ represents the dosage of the ATD that achieves half of the maximum reduction rate (70).

Because patients have varied rates of hyperthyroid progression, we investigate the GIβ dynamics for virtual patients under two different progression rates of hyperthyroidism. We adopt (7-18 pg/ml) and (0.4-4 mU/l) as the normal reference range for the FT4 and TSH respectively (71, 72), and assume U=12.5 pg/mL (the average value of the lower and upper bound of FT4) as the set point of FT4 for the patients we study below. Notably, an apparent deviation from the set point value of FT4 without crossing the bounds of the inter-individual reference range, is sufficient to cause physiological impact on the patients. The American Diabetes Association (ADA) characterizes the fasting glucose levels of euglycemia, pre-diabetes, and diabetes as less than 100 mg/dl, 100-125 mg/dl, and greater than 125 mg/dl, respectively (73). These ranges are contingent on the source and may vary slightly across different labs. In our work, we consider 5 - 20 μU/mL as the reference range for normal fasting insulin and I≥25μ U/mL as the criterion of hyperinsulinemmia (74–76). The model parameters are listed in Table 1 , where the values of h ₂, b ₁, a ₂, b ₂ are chosen to ensure the consistency between simulated hormone dynamics and clinical observations. The remaining parameters are adopted from our previous work (56) where the model was validated by the Pima Inidan data.

We first consider a virtual patient A who is euglycemic and would develop hyperthyroidism with a fast progression rate, as shown in Figure 5B . Figure 5A shows that the increasing amount of excess thyroid hormones drives the gradual and continuous elevation of the glucose level (fasting) before FT4 reaches 19.3 pg/ml. Although the glucose level would be dragged down by 8 mg/dl in 63 days by improved insulin sensitivity, the rising excessive thyroid hormones start deteriorating beta-cells afterwards. The continuous decline of insulin and persistent escalation of glucose follow subsequently. After the gradual elevation of FT4 for 1.7 years, the glucose level of this virtual patient would reach 125mg/dl, the threshold of the diabetic stage. When the FT4 level exceeds its set point value by 124%, a complete beta-cell failure occurs and drives patient A to the late stage of diabetes. The observation of fluctuated blood glucose level is common in clinical settings (77).

Suppose virtual patient B develops hyperthyroidism at a slow rate, spending an extra 17 years more than patient A to reach the FT4 level of 20 pg/ml, as presented in Figure 6B . The time for this patient to develop diabetes would be postponed for 16.7 years compared to patient A, which exhibits the benefit of delayed hyperthyroid progression on slowing the course of diabetes. However, compared to the glucose level of 125 mg/dl that patient A would develop with the FT4 of 24.1 pg/ml, patient B would be in the late stage of diabetes with the same FT4 level. Additionally, the elevated glucose level of patient B is overwhelmingly higher than that of patient A, when their FT4 levels both increase from 20 pg/ml to 24 pg/ml. These results indicate that aside from the amount of excess thyroid hormones, the duration of hyperthyroid exposure affects the severity of hyperglycemia.

We investigate the influence of drug treatment on the GIβ dynamics for the fast progression of hyperthyroidism. The results of timely treated hyperthyroidism are shown in Figure 7 . We assume the ATD drug is administered daily after day 500 with the dose increasing intermittently, that is, the dosage D is designed as a piecewise function, the specific expression of which is described in the caption of Figure 7 . The graphs depict that the initial rise of FT4 level would cause a temporary upsurge of glucose level around the initial point of glucose, which is then quickly counterbalanced by the fast increase of insulin levels. In contrast with the increasing trend of the glucose level in Figure 5A , the administration of 3mg ATD drug after day 500 would result in a sharp reduction of both the glucose level and the insulin level in 10 days.

However, the dose is not enough to counteract the impact of the hyperthyroid factor, and the FT4 level would be elevated again, leading to an upsurge of the glucose level. We then increase the ATD dose to 5 mg for the next 100 days. As a result, the elevated dosage would drive the FT4 down close to its set point value. Subsequently, the glucose level decreases again. Yet this level of treatment is insufficient as the inherent deterioration of the thyroid remains, a condition demonstrated in Figure 5B . To enhance the treatment, an increased dosage is applied every one hundred days until day 800 when the dosage becomes fixed. Consequently, the FT4 level approaches the steady state with 12.9 pg/ml, an improvement over the steady state of Figure 6B . Moreover, the glucose level would be regulated within the normal range after day 1500, a dramatic improvement relative to Figure 5A . In contrast with the massive mortality of beta-cell and the deficient insulin level exhibited in the late period of Figure 5 , the thyroid treatment prevents beta-cell failure and the occurrence of insufficient insulin. Overall, the result highlights that the management of thyroid dysfunction may be of primary consideration for the therapy of diabetic patients.

Notably, our simulation results indicate belatedly treatment of hyperthyroidism fails to reverse diabetes. In Figure 8 , we assume the ATD drug is initiated on day 700 when the FT4 level is 8.2 pg/ml higher than the value on day 500. The dosage regimen in the following days remains the same as in Figure 7 . In this scenario, although the FT4 level can quickly decrease under treatment and maintain a favorable level after day 1000, the ongoing beta-cell failure cannot be halted. This indicates the early diagnosis and timely treatment of hyperthyroidism are key to mitigating diabetes.