A Novel Approach to Numerical Modeling of Metabolic System: Investigation of Chaotic Behavior in Diabetes Mellitus

Although many mathematical models have been presented for glucose and insulin interaction, none of these models can describe diabetes disease completely. In this work, the dynamical behavior of a regulatory system of glucose-insulin incorporating time delay is studied and a new property of the presented model is revealed. This property can describe the diabetes disease better and therefore may help us in deeper understanding of diabetes, interactions between glucose and insulin, and possible cures for this widespread disease.


Introduction
Diabetes, technically called diabetes mellitus, is referred to types of disorders in the metabolic processes of the human body in which the controlling mechanism of sugar level in blood is disrupted. In these cases, insulin, the mainstay controlling element, is either not secreted or its presence is ignored by body cells [1]. There are three types of diabetes: type 1, type 2, and gestational diabetes. Type 1 diabetes is all about insulin. In this type of diabetes, the body's immune system intercepts insulin-releasing cells and destroys them. This type of diabetes accounts for 5 to 10 out of 100 diabetic people. In type 2 diabetes, the body is not able to use insulin in the right way. This type of diabetes is the most common in the world, and around 90 to 95 percent of diabetics have type 2. A third type of diabetes, gestational diabetes, is a temporary condition that occurs during pregnancy. It affects approximately two to four percent of all pregnancies [2][3][4][5].
The function of insulin alters for each organ in the human body, so the effects of environmental factors like stress and nourishing habits may cause blood glucose shift. As observed in various countries, diabetes is discernibly widespread and there is an increasing number of people suffering from this. Hence, the potentially lethal symptoms of the illness necessitate more meticulous treatments and precautionary activities. The essence for such cure procedures is even more accentuated in contemporary hectic life in which people have an increasing penchant to be nourished by artificially cultivated foods and do less exercise. The number of people suffering from this disease was approximately 415 million in 2015 with equal shares of both genders, which accounted for 8.3% of the overall adult population of the world. And nearly 1.5 to 5 million people have died because of diabetes every year between years 2012 and 2015 worldwide [1].
Speaking of the reasons triggering this illness, many elements can cause this irregularity behavior in the body, such as genetic factors inherited through generations that fertilize the body for other factors of the disease to easily disrupt the metabolic system, obesity due to malnutrition and urbanization as consequent of modern lifestyle, side effects of taking specific drugs like glucocorticoids and thyroid hormone, progression of other illnesses, and many other factors which cannot be wholly included [1]. Knowing about the causes of disease enables scientists to develop meditative procedures.
Besides the paramount and distinctive importance of experimental researches for developing effective treatment protocols, studying and developing mathematical models of glucose-insulin bilateral interplay have had an essential role in accelerating the research processes and making breakthroughs in this field by saving both money and time. Conventionally, it was believed that a linear relationship defines the mechanism of glucose-insulin negative feedback system. A linear model for diabetes assumes that the relationship between glucose and insulin concentration could be studied in isolation from other components [6]. In contrast, nonlinear models proposed in previous studies assume that the relationship between components is not always linear [7] and it could depend on initial blood glucose level [8]; moreover, they revealed the fact that statistical properties of the profile in some patients could alter substantially [6,9,10]. In glucose-insulin system, interactions between components are responsible for the overall behavior of the system, which makes this system a complex one. The basic structure of insulin secretion system is a negative feedback controller operating between two elements, namely, the pancreatic β-cells and plasma glucose concentration of the blood contacting these cells. A high level of glucose concentration is acquired, for example, when having a snack which provokes the production and release of insulin leading to a decrease in glucose levels by increasing the consumption rate of the extra sugar or initiation of storage process. On the contrary, if plasma blood is experiencing low levels of glucose concentration, insulin secretion is halted, preventing further declination of blood sugar. In this case, the metabolic system shifts condition from absorptive to postabsorptive [11,12].
As delineated in the preceding paragraph, various mathematical models have been proposed in attempts to simulate the relation between plasma glucose concentration and plasma insulin concentration more accurately, so that scientists will be able to have an elaborate perspective of this metabolic interaction [13][14][15][16][17]. It is wondering that the recently submitted models show some kind of chaotic behavior in the mechanism of malfunctioning metabolic system which is revealed in the current study.
The investigation of chaotic dynamics has attracted the foci of many scientists, and a great deal of effort is put in this field as it has provided a successful method for studying biological systems [18][19][20][21][22][23]. Moreover, this novel vantage point of studying biological phenomena has made revolutionary effects on developing biological system models [24][25][26][27].
Because of the complexity of the system, the model that has been studied in this paper is a nonlinear model. The nonlinear model that we study reflects the relationship between injected insulin and blood glucose response. The studies about variation in the blood glucose indicate a chaotic component.
In the second section, the dynamical properties of the last presented model for glucose and insulin concentration are investigated. Eventually, conclusion remarks are given in Section 4.

Mathematical Model
In 1964, Ackerman et al. [16] proposed a linear model for glucose tolerance test consisting of two ordinary differentials (1), as demonstrated below.
1 where x t is the insulin concentration and y t is the blood glucose concentration. a 1 y t is the rate of increase in insulin concentration due to increase in glucose concentration, a 2 x t represents the rate of insulin reduction, a 3 y t represents the rate of glucose reduction independent to insulin, and a 4 x t represents the rate of glucose removal due to insulin secretion. C 1 and C 2 are positive constants, and I t is the rate of increase in blood glucose concentration due to absorption in the gastrointestinal system. In 1987, Bajaj et al. [15] proposed a nonlinear mathematical model for glucose-insulin feedback system which incorporated β-cell kinetics. The mathematical relationships for the model are formulated as shown in where x t and y t represent the insulin and glucose concentrations, respectively, and z t represents the number of β-cells. It has been discovered that β-cells have an essential role in regulating glucose and insulin concentration. Recent studies indicate two delays in glucose-insulin feedback control system [28][29][30][31][32]. Two important time lags can be noticed in the system, the lag in insulin secretion in response to an increase in blood glucose concentration, τ g , and hepatic glucose response lag, τ i . In the current research, we study a nonlinear mathematical model for glucose-insulin feedback control system by incorporating the enhanced delay differential equations embracing β-cells proposed by the model presented by Sarika et al. [28]. The modified model is the compound of the model proposed by Bajaj et al. [15] and the one suggested by Sarika et al. [28]. The resulted model is the one presented by 2 Complexity Chuedoung et al. [11] (Figure 1), which can be presented in delay differential equations as follows: where x t is the insulin concentration, y t is the glucose concentration, z t is the number of β-cells, andŷ is the difference between glucose fasting level and its basal level. τ g is the delay in insulin secretion in response to blood glucose level increase based on clinical evidence reported by Palumbo et al. [33], and τ i is the delay in glucose drop due to increased insulin level based on clinical evidence reported by Prager et al. [34]. r 1 y t − τ g z t − τ g shows the increase in insulin concentration in response to blood glucose increase with the time delay τ g . r 2 x is the rate of insulin decrease independent of glucose, and c 1 z t − τ g is the increase of insulin level secreted by β-cells and is independent from other components. System (3) considers two time lags in insulin-glucose regulatory system; therefore, it is more realistic and is capable of showing the behavior of insulin-glucose regulatory system in different time delays. Previous models cannot display the behavior of aforementioned biological system with respect to time delays. According to the model presented by Molnar et al. [17], if insulin secretion decreases to 1/N of the number of β-cells, designated by N, due to a reduction, then the blood glucose increases until insulin levels are restored to nearly normal standards. So the blood glucose level is a function of the β-cells' capacity N/n. N is the normal number of β-cells. R 4 x t − τ i is the rate of glucose reduction in response to insulin secretion with the time delay τ i . T is the total density of β-cells, and the term R 5 y −ŷ T − z represents the increase in dividing β-cells caused by the interaction between blood glucose above the fasting level and the nondividing β-cells. The term R 6 z T − z represents the increase in z due to interaction between dividing and nondividing β-cells, and the term R 7 z represents the reduction in z due to its current level.

Results and Discussion
Based on the study by Chuedoung et al. [11], the mentioned model shows the different behaviors for different parameters. The proposed model comprises a number of parameters that their values are essential in changing the behavior of the system.
In current research, the new capability of the mentioned model is revealed. By increasing the insulin secretion delay by β-cells (τ g ), the system behaves in a chaotic way ( Figure 2). Figure 2 is the bifurcation diagram of the system for the different values of τ g . Figure 2 shows that if there is more delay on insulin secretion, insulin cannot track glucose and the concentration of blood glucose rises which results in diabetes disorder. Technically speaking, the time lag of insulin response in glucose-insulin negative feedback controlling mechanism is shown to be the main reason for this disease. Computer simulation of  [11].
With the increase in glucose response time caused by insulin secretion (τ i ), the system behaves chaotically (Figure 3). Figure 3 is the bifurcation diagram of the system for different values of τ i . It shows that if there is more delay on response of glucose to insulin secretion, the system behaves in a chaotic manner. Note that the goal of this study is not to investigate the parameters quantitatively, but rather to show that a minute change in the quantities of parameters of the model can bring about changes in the behavior of the system. Computer simulation of Figure 3 is done by the following parameters r 1 = 0 472, r 2 = 0 25, R 3 = 0 82, R 4 = 0 6, R 5 = 0 3, R 6 = 0 3, R 7 = 0 2, c 1 = 0 1, C 2 = 0 8,ŷ = 1 42, T = 1 5, N = 1 27, and τ g = 0 56. The values of the parameters are the same as the values presented by Chuedoung et al. [11].
In the present study, we used the insulin-glucose model involving β-cells presented by Chuedoung et al. [11] and the effect of delays on insulin-glucose model have been investigated. The system is stable for small delays, and when the delays increase, the system exhibits chaotic behavior. According to the claim made by Bertram and Pernarowski [35], 1-2 min lag, representative of insulin secretion, is a common incident after bath application of glucose in islet electrical activity when investigating islet porosity and the permeability of a surrounding layer of acinar cells on the time required for glucose to diffuse through an isolated pancreatic islet of Langerhans and reach an equilibrium. And, based on a report by Forrest et al. [36], instantaneous insulin reflection was recorded in 14 out of the 20 monitored Jamaican children rehabilitated from malnutrition. The response time was about 1 minute for them; whereas, this delay ranged from 5 to 10 minutes for the other children. Hence, these observations support our claim about the range of τ g and τ i .

Adaptive Sliding Mode Control
A new sliding mode control scheme for a class of uncertain time-delay chaotic systems is proposed in [37]. It is shown that a linear time-invariant system with the desired system dynamics is used as a reference model for the output of a time-delay chaotic system to track. Chaos control for scalar-delayed chaotic systems using sliding mode control strategy is achieved in [38]. Sliding surface design is based on delayed feedback controller, and it is shown that the proposed controller can achieve stability for an arbitrary unstable fixed point (UPF) or unstable periodic orbit (UPO) with an arbitrary period.
In this section, we design the adaptive sliding mode controllers to suppress the chaotic oscillations in the model presented in (3). For the uncertainties, we assume that the parameters r 1 , r 2 , c 1 , and c 2 are unknown. The entire control algorithm is designed with the delay elements as described in (3), and hence the sliding surface initialization is accounted with the respective time delays. Let us redefine the model in (3) with the controllers u i , where i = x, y, z as given in We define the integral sliding mode surface as x τ dτ, s y = y + k y t 0 y τ dτ, The sliding surface dynamics can be derived as s y = y + k y y, The parameter estimation errors are defined as e r 1 =r 1 − r 1 , e r 2 =r 2 − r 2 , The first derivatives of the estimation errors are Consider the following Lyapunov function: The first derivative of the Lyapunov candidate function is V = s 1 s 1 + s 2 s 2 + s 3 s 3 + e r 1 e r 1 + e r 2 e r 2 + e c 1 e c 1 + e c 2 e c 2 10 Applying (4), (6), and (8) in (10), we have By introducing uncertainties without changing the definition in (11), x t − τ i + c 2 + c 2 − c 2 + u y + k y y + s z R 5 y − y T − z + R 6 z T − z − R 7 z + u z + k z z + e r 1 r 1 + e r 2 r 2 + e c 1 c 1 + e c 2 c 2 12

Complexity
After some mathematical simplifications, let us define the adaptive sliding mode controllers as The parameter estimate laws can be defined as Using (13) and (14) in (12), we simplify the Lyapunov candidate function dynamics to As ρ i and η i are positive for i = x, y, z, the Lyapunov first derivative (15) is a negative definite function which infers that the controller is stable as per the theorem discussed in [39,40] and is valid for any bounded initial conditions. For numerical simulations, the initial conditions are taken as 2.6, 1, and 0.825 and the sliding surface initial conditions are defined as −2.6, −1, and −0.825 with the time delays  Figure 6 shows the estimated parameters with parameter update laws and controllers in action at t = 0 s.

FPGA Implementation
Implementation of chaotic and hyperchaotic systems using Field Programmable Gate Arrays (FPGA) has been widely investigated [41][42][43]. Chaotic random number generators have been implemented in FPGA for applications in image cryptography [44]. FPGA-implemented Duffing oscillatorbased signal detectors have been proposed by Rashtchi et al. [45]. Digital implementations of chaotic multiscroll 6 Complexity attractors have been extensively investigated [41,46]. Memristor-based chaotic system and its FPGA circuits have been proposed by Ya-Ming et al. [47]. A FPGA implementation of fractional order chaotic system using approximation method has been investigated by Rajagopal et al. [48][49][50].
In this section, we implement a circuit for the model (3) by FPGA. To the best of our knowledge, only a few literatures [51,52] have implemented delay chaotic systems. However, those works discuss about indirect realizations which will increase the time slack factor as the programs run sequentially on the processor. But we use a direct realization, and hence the power utilization and the time slack delays are reduced. For the design of delay chaotic model (3), first, we configure the available built-in blocks of the System generator toolbox. The Add/Sub blocks are configured with zero latency and 32/16 bit fixed point settings. The delays are introduced by an additional delay block introduced with the defined time delays as in (3). The output of the block is configured to rounded quantization in order to reduce the bit latency. Then, we design the integer order integrator which is not a readily available block in the System Generator. Hence, we implement the integrators using the mathematical relation dx i /dt = lim h→0 x i n + 1 − x i n /h and the value of h is taken as 0 001.

Complexity
The initial conditions are fed into the forward register. Figure 7 shows the overall RTL schematics of the system with time delays. Figure 8 shows the RTL schematics of the implemented delay (3). Figure 9 shows the Xilinx RTL schematics of the controllers with sliding surfaces and parameter estimate laws. Figures 10 and 11 show the controlled states of the delay systems and estimated parameters with the control in action at t = 0 5 s, respectively.

Discussion and Conclusion
By studying presented mathematical models for glucoseinsulin interaction, according to the value of parameters, a chaotic model for describing the glucose-insulin regulatory system was found. In the present study, it is expected to observe periodic behavior in the proposed system under normal metabolic conditions and chaotic behavior under abnormal metabolic conditions. It is noteworthy to say that the chaotic behavior of a system is a sign of a faulty condition in the biological systems [18][19][20][21][22][23].
The effect of two time delays on glucose-insulin regulatory system was investigated. Two main results of this study are listed as below.
(i) If the time lag of insulin response to glucose increases, system exits from periodic region and enters to chaotic region, and if there is more delay on response of glucose to insulin secretion, the system displays chaotic manner, which was in line with previous studies [35,36]  8 Complexity (ii) If there is more delay on the response of glucose to insulin secretion, the behavior of system alters from periodic to chaotic The proposed system can explain the interaction between glucose and insulin concentration in both normal and abnormal (diabetes disease) situations. Also, a control method was investigated with the hope of possible clinical applications.

Conflicts of Interest
The authors declare that they have no conflicts of interest.